!Last change: CYCO 28 Dec 104 2:49 pm
!modules are highlighted by ***, Subroutines by ===
!****************************
module NumKind
!****************************
!This module defines the kind of integer and real number
!Every module, subroutine and function must use this module.
implicit none
integer(kind(1)),parameter::ikind=kind(1),rkind=kind(0.D0)
real(rkind),parameter::Zero=0.D0, One=1.D0, Two=2.D0, Three=3.D0,&
Four=4.D0, Five=5.D0, Six=6.D0, Seven=7.D0,&
Eight=8.D0, Night=9.D0, Ten=10.D0, Twelve=12.D0
!Integer from 0-10 and 12 relating to any calculation should be written as above,
!for the convenience of changing the precision of the analysis.
end module NumKind
!****************************
module TypeDef
!****************************
use NumKind
implicit none
!Every joint should be stored as the following format:
type::typ_Joint
real(rkind) :: x,y !The two coordicates of a joint
integer(ikind) :: GDOF(3) !The three Global DOF of a joint
end type typ_Joint
!Every Element should be stored as the following format:
type::typ_Element
!The index number of the starting and the end joint of an element
integer(ikind) ::JointNo(2)
!The Element Location Vector
integer(ikind) ::GlbDOF(6)
!The Element properties
real(rkind) ::Length,CosA,SinA,EI,EA
end type typ_Element
!Every Joint Load should be stored as the following format:
type::typ_JointLoad
!Here load can be force, moment, or displacement.
integer(ikind) :: JointNo !The no. of the Joint where the Load is applied
integer(ikind) :: LodDOF !The direction of the load
real(rkind) :: LodVal !The magnitude of the load
end type typ_JointLoad
!Every Element Load should be stored as the following format:
type::typ_Elemload
!Here load can be force, pressure, moment or temperature changes
integer(ikind) :: ElemNO !the no. of the element where the Load is applied
integer(ikind) :: Indx !The type index of the Load
real(rkind) :: pos,LodVal !The position and magnitude of the load
!Here the Pos refers to the relative distance of the length of the element.
end type typ_Elemload
contains
!========================
subroutine SetElemProp(Elem,Joint)
!========================
!Calculate the basic properties of an element according to the
!joints' information of the element
type(typ_Element),intent(in out) ::Elem(:)
type(typ_Joint),intent(in) ::Joint(:)
integer(ikind) ::ie,Nelem
real(rkind) ::x1,x2,y1,y2
Nelem=size(Elem,1)
do ie=1,Nelem
x1 = Joint(Elem(ie)%JointNo(1))%x
x2 = Joint(Elem(ie)%JointNo(2))%x
y1 = Joint(Elem(ie)%JointNo(1))%y
y2 = Joint(Elem(ie)%JointNo(2))%y
Elem(ie)%Length = sqrt((x2-x1)**Two+(y2-y1)**Two)
Elem(ie)%CosA = (x2-x1)/Elem(ie)%Length
Elem(ie)%SinA = (y2-y1)/Elem(ie)%Length
Elem(ie)%GlbDOF(1:3) = Joint(Elem(ie)%JointNo(1))%GDOF(1:3)
Elem(ie)%GlbDOF(4:6) = Joint(Elem(ie)%JointNo(2))%GDOF(1:3)
end do
return
end subroutine SetElemProp
!========================
subroutine TransMatrix(ET,CosA,SinA)
!========================
!Creat the transform matric of the element ET
real(rkind),intent(out):: ET(:,:)
real(rkind),intent(in) :: CosA,SinA
! ET could be 2x2, 3x3 or 6x6 depending on size(ET)
ET=Zero
ET(1,1:2)=(/CosA,SinA/)
ET(2,1:2)=(/-SinA,CosA/)
if(size(ET,1)>2) ET(3,3)=One
if(size(ET,1)>3) ET(4:6,4:6)=ET(1:3,1:3)
return
end subroutine TransMatrix
end module TypeDef
!************************
module BandMat
!************************
use NumKind
use TypeDef,only:typ_Element
implicit none
private
public::SetMatBand,DelMatBand,VarBandSolv
type,public::typ_Kcol
real(rkind),pointer::row(:)
end type typ_Kcol
contains
!========================
subroutine SetMatBand(Kcol,Elem)
!========================
type(typ_Kcol),intent(in out):: Kcol(:)
type(typ_Element), intent(in):: Elem(:)
integer(ikind) :: minDOF,ELocVec(6)
integer(ikind) :: ie,j,NGlbDOF,NElem
integer(ikind) :: row1(size(Kcol,1))
NGlbDOF = size(Kcol,1)
NElem = size(Elem,1)
row1 = NGlbDOF
!?·?¨?÷??????????·?????×é row1(:)??
do ie=1,Nelem
ELocVec=Elem(ie)%GlbDOF
minDOF =minval(ELocVec,mask=ELocVec>0)
where(ELocVec>0)
row1(ELocVec)=min(row1(ELocVec),minDOF)
end where
end do
!???÷????°????í·???????????????
do j=1,NGlbDOF
allocate(Kcol(j)%row(row1(j):j))
Kcol(j)%row=Zero !????
end do
return
end subroutine SetMatBand
!========================
subroutine DelMatBand(Kcol)
!========================
type(typ_Kcol),intent(in out):: Kcol(:)
integer(ikind) :: j,NGlbDOF
NGlbDOF = size(Kcol,1)
do j=1,NGlbDOF
deallocate (Kcol(j)%row)
end do
return
end subroutine DelMatBand
!========================
subroutine VarBandSolv(Disp,Kcol,GLoad)
!========================
type (typ_KCol),intent(in out) :: Kcol(:)
real (rkind),intent(out) :: Disp(:)
real (rkind),intent(in) :: GLoad(:)
integer (ikind) :: i,j,row1j,row_1,NCol
real(rkind) :: s,Diag(size(Kcol,dim=1))
NCol=size(Kcol,1)
Diag(:)=(/(Kcol(j)%row(j),j=1,NCol)/)
do j=2,NCol
row1j=lbound(Kcol(j)%row,1)
do i=row1j,j-1
row_1=max(row1j,lbound(Kcol(i)%row,1))
s=sum(Diag(row_1:i-1)*Kcol(i)%row(row_1:i-1)*Kcol(j)%row(row_1:i-1))
Kcol(j)%row(i)=(Kcol(j)%row(i)-s)/Diag(i)
end do
s=sum(Diag(row1j:j-1)*Kcol(j)%row(row1j:j-1)**Two)
Diag(j)=Diag(j)-s
end do
Disp(:) = GLoad(:)
!...[ 6-5-3 ?????ú???????? GP ???? Disp ]
do j=2,NCol
row1j=lbound(Kcol(j)%row,1)
Disp(j)=Disp(j)-sum(Kcol(j)%row(row1j:j-1)*Disp(row1j:j-1))
end do
!...[ 6-5-4 ?????ú???????? GP ???? Disp ]
Disp(:)=Disp(:)/Diag(:)
do j=NCol,1,-1
row1j=lbound(Kcol(j)%row,1)
Disp(row1j:j-1)=Disp(row1j:j-1)-Disp(j)*Kcol(j)%row(row1j:j-1)
end do
return
end subroutine VarBandSolv
end module BandMat
!****************************
module DispMethod
!****************************
!The displacement method in structure machemism
use NumKind
use TypeDef
use BandMat
implicit none
contains
!========================
subroutine SolveDisp(Disp,Elem,Joint,JLoad,ELoad)
!========================
real(rkind),intent(out) :: Disp(:)
type(typ_Element),intent(in) :: Elem(:)
type(typ_Joint) ,intent(in) :: Joint(:)
type(typ_JointLoad),intent(in) :: JLoad(:)
type(typ_Elemload),intent(in) :: ELoad(:)
type(typ_Kcol),allocatable :: Kcol(:)
real(rkind),allocatable :: GLoad(:)
integer(ikind) :: NGlbDOF
NGlbDOF=size(Disp,dim=1)
allocate(GLoad(NGlbDOF))
allocate(Kcol(NGlbDOF))
call SetMatBand(Kcol,Elem)
call GLoadVec(GLoad,Elem,JLoad,ELoad,Joint)
call GStifMat(Kcol,Elem)
call VarBandSolv(Disp,Kcol,GLoad)
call DelMatBand(Kcol)!free the memory space occupied by Kcol
deallocate(GLoad) !free the memory space occupied by GLoad
return
end subroutine SolveDisp
!========================
subroutine GStifMat(Kcol,Elem)
!========================
!Creat the Globle Stiffness Matrix for the problem
type(typ_Kcol),intent(in out) :: Kcol(:)
type(typ_Element),intent(in) :: Elem(:)
integer (ikind) :: ie,j,JGDOF,NElem
real (rkind) :: EK(6,6),ET(6,6)
integer (ikind) :: ELocVec(6)
NElem=size(Elem,1)
do ie=1,NElem
!Create the Element Stiffness Matrix
call EStifMat(EK,Elem(ie)%Length,Elem(ie)%EI,Elem(ie)%EA)
!Create the Transform Matrix
call TransMatrix(ET,Elem(ie)%CosA,Elem(ie)%SinA)
!Transform the Element stiffness Matrix to Globle coordinate system
EK=matmul(transpose(ET),matmul(EK,ET))
!Get the Element Location Vector
ELocVec=Elem(ie)%GlbDOF
!Integrate the element stiffness matrix to the Global Stiffness matrix
!according to their Element Location Vectors.
do j=1,6
JGDOF=ELocVec(j)
if(JGDOF==0) cycle
where (ELocVec<=JGDOF.and.ELocVec>0)
Kcol(JGDOF)%row(ELocVec)=Kcol(JGDOF)%row(ELocVec)+EK(:,j)
end where
end do
end do
return
end subroutine GStifMat
!========================
subroutine GLoadVec(GLoad,Elem,JLoad,ELoad,Joint)
!========================
real(rkind), intent(out):: GLoad(:)
type(typ_Element), intent(in) :: Elem(:)
type(typ_JointLoad),intent(in) :: JLoad(:)
type(typ_Elemload), intent(in) :: ELoad(:)
type(typ_Joint), intent(in) :: Joint(:)
integer(ikind) :: i !temporary Loop variable
integer(ikind) :: NJLoad !Number of Joint Load
integer(ikind) :: NELoad !Number of Element Load
integer(ikind) :: ELocVec(6),ve(6),GDOF(3),vj,LodDOF
real(rkind) :: F0(6) !Fix End Forces
real(rkind) :: ET(6,6)
integer(ikind) :: ElemNo,JointNo
NJLoad = size(JLoad,1)
NELoad = size(ELoad,1)
GLoad = Zero !Initialize
ve(:) = (/1,2,3,4,5,6/) !An assistant vector
!Deal with every element load
do i=1,NELoad
ElemNo=ELoad(i)%ElemNo
ELocVec(:)=Elem(ElemNo)%GlbDOF(:) !Get the Element Location Vector
F0=Zero !Initialize
!Calculate the element fixend forces vector F0(6) under the local coordinate system
call EFixendF(F0,ELoad(i)%Indx,ELoad(i)%pos,ELoad(i)%LodVal,Elem(ElemNo))
!Create the Transform Matrix
call TransMatrix(ET,Elem(ElemNo)%CosA,Elem(ElemNo)%SinA)
!Transform the fixend forces vector to Globle coordinate system
F0=matmul(transpose(ET),F0)
!Integrate F0 to the Global Load Vector
where (ELocVec>0)
GLoad(ELocVec)=GLoad(ELocVec)-F0(ve)
end where
end do
!Deal with every joint load
do i=1,NJLoad
JointNo = JLoad(i)%JointNo !The joint where the load is applied
GDOF(:) = Joint(JointNo)%GDOF(:) !The Global DOF of the joint
LodDOF = JLoad(i)%LodDOF !The direction of the load (local coordinate system)
vj=GDOF(LodDOF) !The direction of the load (global coordinate system)
!Integrate JLoad to the Global Load Vector
GLoad(vj)=GLoad(vj)+JLoad(i)%LodVal
end do
return
end subroutine GLoadVec
!========================
subroutine EStifMat(EK,Elen,EI,EA)
!========================
!Create the element Stiffness Matrix EK
real(rkind),intent(out) :: EK(6,6)
real(rkind),intent(in) :: Elen,EI,EA
integer(ikind) :: i,j
EK=Zero !Initialize
!Generate the upper triangle part of EK
EK(1,1) = EA/Elen
EK(1,4) = -EA/Elen
EK(2,2) = Twelve*EI/Elen**Three
EK(2,3) = Six*EI/Elen**Two
EK(2,5) = -Twelve*EI/Elen**Three
EK(2,6) = Six*EI/Elen**Two
EK(3,3) = Four*EI/Elen
EK(3,5) = -Six*EI/Elen**Two
EK(3,6) = Two*EI/Elen
EK(4,4) = EA/Elen
Ek(5,5) = Twelve*EI/Elen**Three
EK(5,6) = -Six*EI/Elen**Two
EK(6,6) = Four*EI/Elen
!Mirror the upper triangle part to get the complete matrix
do j=1,6
do i=j+1,6
EK(i,j)=EK(j,i)
end do
end do
return
end subroutine EStifMat
!========================
subroutine EFixendF(F0,Indx,a,q,Elem)
!========================
!Calculate the fixend forces vector F0 of an element
real(rkind), intent(out):: F0(6)
integer(ikind),intent(in) :: Indx
real(rkind), intent(in) :: a !the relative distance
real(rkind), intent(in) :: q !the magnitude of the load
type(typ_Element),intent(in) :: Elem
real ::l
l=Elem%length
select case (Indx)
case (1) !Transverse uniform pressure
F0(1) = Zero
F0(2) = -q*(a*l)/Two*(Two-Two*a**Two+a**Three)
F0(3) = -q*(a*l)**Two/Twelve*(Six-Eight*a+Three*a**Two)
F0(4) = Zero
F0(5) = -q*(a*l)*a**Two/Two*(Two-a)
F0(6) = q*(a*l)**Two*a/Twelve*(Four-Three*a)
case (2) !Transverse concentrated force
F0(1) = Zero
F0(2) = -q*(One-Two*a+a**Two)*(One+Two*a)
F0(3) = -q*(a*l)*(One-Two*a+a**Two)
F0(4) = Zero
F0(5) = -q*a**Two*(Three-Two*a)
F0(6) = q*a**Two*(One-a)*l
case (3) !axial displacement of the supports
if(a<One/Two)then
F0(1)=Elem%EA*q/l
F0(4)=-F0(1)
else
F0(1)=-Elem%EA*q/l
F0(4)=-F0(1)
end if
F0(2) = Zero
F0(3) = Zero
F0(5) = Zero
F0(6) = Zero
case (4) !transvers displacement of the supports
if(a<One/Two)then
F0(2) = Twelve*Elem%EI*q/l**Three
F0(3) = Six*Elem%EI*q/l**Two
F0(5) = -F0(2)
F0(6) = F0(3)
else
F0(2) = -Twelve*Elem%EI*q/l**Three
F0(3) = -Six*Elem%EI*q/l**Two
F0(5) = -F0(2)
F0(6) = F0(3)
end if
F0(1) = Zero
F0(4) = Zero
end select
return
end subroutine EFixendF
!========================
subroutine ElemDisp(EDisp,ie,Disp,Elem)
!========================
!Get the Element Displacement from the global result Disp()
!in the Global Coordinate System
real(rkind), intent(out) :: EDisp(:) !Element displacement to be generated
real(rkind), intent(in) :: Disp(:) !Global displacement solve
integer(ikind), intent(in) :: ie !Index of the element to be dealt with
type(typ_Element),intent(in) :: Elem(:) !All the Elements
integer(ikind) :: i !loop variable
do i=1,6
if(Elem(ie)%GlbDOF(i)==0)then
EDisp(i)=Zero
else
EDisp(i)=Disp(Elem(ie)%GlbDOF(i))
end if
end do
return
end subroutine ElemDisp
!========================
subroutine ElemForce(EForce,ie,Disp,Elem,ELoad)
!========================
!Get the Element end forces from the global result
!in the Local Coordinate System
!Here we have to repeat the function contained in the subroutine ElemDisp
!in consideration of independence between the calculations of displacement and forces
real(rkind),intent(out):: EForce(:)
real(rkind),intent(in) :: Disp(:)
integer(ikind) ,intent(in) :: ie
type(typ_Element) ,intent(in) :: Elem(:)
type(typ_ElemLoad),intent(in) :: Eload(:)
real(rkind) :: EDisp(6)
real(rkind) :: ET(6,6),EK(6,6)
real(rkind) :: EP(6) !Fix end Forces
integer(ikind) :: NELoad !Number of Element Load
integer(ikind) :: i !loop variable
integer(ikind) :: ElemNo !Where the element load is applied
!Creat a transform matrix
call TransMatrix(ET,ELem(ie)%CosA,Elem(ie)%SinA)
!Get the element displacement
do i=1,6
if(Elem(ie)%GlbDOF(i)==0)then
EDisp(i)=Zero
else
EDisp(i)=Disp(Elem(ie)%GlbDOF(i))
end if
end do
!Creat the element stiffness matrix
call EStifMat(EK,Elem(ie)%Length,Elem(ie)%EI,Elem(ie)%EA)
!Calculate the element end forces caused by displacement
EForce=matmul(EK,matmul(ET,EDisp))
!Calculate the element end forces caused by element load
!Which is the same as end fix forces
NELoad = size(ELoad,1)
EP = Zero !Initialize
do i=1,NELoad
if(ELoad(i)%ElemNo==ie)then
ElemNo=ELoad(i)%ElemNo
!Calculate the element fixend forces vector F0(6) under the local coordinate system
call EFixendF(EP,ELoad(i)%Indx,ELoad(i)%pos,ELoad(i)%LodVal,Elem(ElemNo))
!Add the two part to get the whole element end forces
EForce(:) = EForce(:)+EP(:)
end if
end do
return
end subroutine ElemForce
end module DispMethod
!*****************************
program SM_90 !main program
!*****************************
use NumKind
use TypeDef
use DispMethod
implicit none
type(typ_Element) ,allocatable :: Elem(:)
type(typ_Joint) ,allocatable :: Joint(:)
type(typ_JointLoad),allocatable :: JLoad(:)
type(typ_Elemload) ,allocatable :: ELoad(:)
real(rkind) ,allocatable :: Disp(:)
integer(ikind) :: ProbType
integer(ikind) :: NElem,NJoint,NGlbDOF,NJLoad,NELoad
integer(ikind) :: i,ie
call InputData() !See below
if(ProbType==1)then
call SetElemProp(Elem,Joint)
call SolveDisp(Disp,Elem,Joint,JLoad,Eload)
end if
call OutputResutls() !See below
stop
contains
!========================
subroutine InputData()
!========================
open(5,file="SMCAI90.IPT",action="READ",position="REWIND",status="OLD")
read(5,*) ProbType
read(5,*) NElem,NJoint,NGlbDOF,NJLoad,NELoad
allocate(Elem(NElem))
allocate(Joint(NJoint))
allocate(JLoad(NJLoad))
allocate(ELoad(NELoad))
allocate(Disp(NGlbDOF))
Disp=Zero
read(5,*) (Joint(i),i=1,NJoint)
read(5,*) (Elem(ie)%JointNo,Elem(ie)%EA,Elem(ie)%EI,ie=1,NElem)
if(NJLoad>0)read(5,*)(JLoad(i),i=1,NJLoad)
if(NELoad>0)read(5,*)(ELoad(i),i=1,NELoad)
return
end subroutine InputData
!========================
subroutine OutputResutls()
!========================
real(rkind) :: EDisp(6),EForce(6),EDispLoad(6)
real(rkind) :: ET(6,6)
integer(ikind) :: ie
open(8,file="SMCAI90.OUT",action="WRITE",position="REWIND",status="REPLACE")
write(8,*) "10 0"
do ie=1,size(Elem)
call ElemDisp(EDisp,ie,Disp,Elem)
EDispLoad = Zero !Displacement Load Initialize
!Create the Transform Matrix
call TransMatrix(ET,Elem(ie)%CosA,Elem(ie)%SinA)
!get the EDispLoad(6) under the local coordinate system
do i=1,NELoad
if(ELoad(i)%ElemNo==ie)then
if(ELoad(i)%Indx==3)then
if(ELoad(i)%pos==0)then
EDispLoad(1)=EDispLoad(1)+ELoad(i)%LodVal
else
EDispLoad(4)=EDispLoad(4)+ELoad(i)%LodVal
end if
else if(ELoad(i)%Indx==4)then
if(ELoad(i)%pos==0)then
EDispLoad(2)=EDispLoad(2)+ELoad(i)%LodVal
else
EDispLoad(5)=EDispLoad(5)+ELoad(i)%LodVal
end if
end if
end if
end do
!Transform the EDispLoad(6) to global coordinate system
EDispLoad=matmul(transpose(ET),EDispLoad)
!Add the two part of displacement together
EDisp=EDisp+EDispLoad
write(8,*) EDisp
end do
do ie=1,size(Elem)
call ElemForce(EForce,ie,Disp,Elem,ELoad)
write(8,*) -EForce(1),EForce(2:6)
end do
deallocate(JLoad)
deallocate(Joint)
deallocate(Elem)
deallocate(ELoad)
deallocate(Disp)
return
end subroutine OutputResutls
end program SM_90
转载自原文链接, 如需删除请联系管理员。
原文链接：结构力学求解器v2.0 源码奉献 可以运行，转载请注明来源！