!===========================================================
! 非线性分析程序设计范例
! 利用Newton-Raphson法求解平面弹塑性桁架
！程序说明
! 程序作者：陆新征
! 北京清华大学土木工程系
! 2006.5.
!===========================================================

module TypDef
implicit none
integer (kind(1)),parameter ::ikind=(kind(1))
integer (kind(1)),parameter ::rkind=(kind(0.D0))

type :: Typ_GValue ! 整体控制参数
integer (ikind) :: NNode, NElem, NLoad, NSupport ! 节点，单元，荷载，支座数目
integer (ikind) :: NStep, CurrentStep ! 总加载步数, 当前加载步数
real (rkind) :: Tol,Error ! 误差容限, 当前误差
end type Typ_GValue

type :: Typ_Node ! 定义节点
real(rkind) :: Coord(2) ! 节点坐标
real(rkind) :: Disp(2),dDisp(2) ! 节点总位移，节点位移增量
end type Typ_Node

type :: Typ_Truss ! 定义桁架截面
integer(ikind) :: NodeNo(2) ! 节点编号
real(rkind) :: E, A, L ! 弹性模量，截面积，长度
real(rkind) :: fy, fy_temp, Eh ! 屈服强度，，屈服强度临时变量，强化模量
real(rkind) :: Et ! 迭代刚度
real(rkind) :: T(4,4) ! 坐标转化矩阵
real(rkind) :: B(1,2) ! 几何矩阵
real(rkind) :: K(4,4) ! 单元刚度矩阵
real(rkind) :: EPS, dEPS, SIG ! 迭代开始时应变，应变增量，应力
real(rkind) :: STRESS ! 迭代结束时应力
real(rkind) :: R(4) ! 单元节点反力
end type Typ_Truss

type :: Typ_Load ! 定义荷载
integer (ikind) :: NodeNo ! 作用节点
real(rkind) :: Value(2) ! 荷载大小
end type Typ_Load

type :: Typ_Support ! 定义支座
integer (ikind) :: NodeNo ! 作用节点
integer (ikind) :: DOF ! 所约束的自由度
end type Typ_Support

contains

subroutine Get_dDisp(GValue, Node, GDisp) ! 计算节点位移增量
type(typ_GValue) :: GValue
type(Typ_Node) :: Node(:)
real(rkind) :: GDisp(:)
integer (ikind) :: i
do i=1, size(Node)
Node(i)%dDisp(1:2)=Node(i)%dDisp(1:2)+GDisp(i*2-1:i*2)
end do
return
end subroutine Get_dDisp

subroutine Get_Stress(Elem) ! 计算应力增量和迭代弹性模量
type(Typ_Truss) :: Elem
real(rkind) :: S_trial, dE_pla, x
Elem%Et=Elem%E
S_trial=Elem%SIG+Elem%dEPS*Elem%E ! 得到弹性应力增量
if(abs(S_trial)>Elem%fy) then ! 如果超过屈服应力
dE_pla=(abs(S_trial)-Elem%fy)*sign(1.d0,S_trial)/Elem%E ! 得到塑性应变增量
Elem%Stress=Elem%fy*sign(1.d0,S_trial)+dE_pla*Elem%Eh
Elem%fy_temp=abs(Elem%Stress) ! 更新屈服应力
if(Elem%dEPS.ne.0.d0) then
Elem%Et=Elem%Eh
end if
else
Elem%Stress=S_trial
end if
return
end subroutine Get_Stress

subroutine Get_EK(Elem) ! 组装单元刚度矩阵
type(Typ_Truss) :: Elem
real(rkind) :: K(2,2)
K=matmul(transpose(Elem%B),Elem%B)*Elem%Et*Elem%L*Elem%A
Elem%K=0.d0
Elem%K(1,1)=K(1,1);
Elem%K(3,3)=K(2,2);
Elem%K(1,3)=K(1,2);
Elem%K(3,1)=K(2,1)
Elem%K=matmul(transpose(Elem%T),matmul(Elem%K,Elem%T)) ! K=T'KT
return
end subroutine Get_EK

subroutine Get_R(Elem) ! 得到单元反力
type(Typ_Truss) :: Elem
real(rkind) :: R0(2,1), R1(4) ! 单元坐标系下的反力
R0=transpose(Elem%B)*Elem%Stress*Elem%L*Elem%A ! R=B*SIG*L*A
R1=0.d0;
R1(1)=R0(1,1); R1(3)=R0(2,1);
Elem%R=matmul(transpose(Elem%T),R1) ! R=T'R'
return
end subroutine Get_R

subroutine Get_EPS(Node,Elem) ! 得到应变增量
type(typ_Node) :: Node(:)
type(typ_Truss) :: Elem
real(rkind) :: D(4)
D(1:2)=Node(Elem%NodeNo(1))%dDisp(1:2)
D(3:4)=Node(Elem%NodeNo(2))%dDisp(1:2)
D=matmul(Elem%T,D) ! D'=TD
Elem%dEPS=D(1)*Elem%B(1,1)+D(3)*Elem%B(1,2)
return
end subroutine Get_EPS

subroutine Intial_Elem(Node,Elem) ! 初始化单元，得到单元坐标变化矩阵和几何矩阵
type(typ_Node) :: Node(:)
type(typ_Truss) :: Elem(:)
integer(ikind) :: ElemNo
real(rkind) :: x1,x2,y1,y2,cosA,sinA
do ElemNo=1,size(Elem)
x1=Node(Elem(ElemNo)%NodeNo(1))%Coord(1)
x2=Node(Elem(ElemNo)%NodeNo(2))%Coord(1)
y1=Node(Elem(ElemNo)%NodeNo(1))%Coord(2)
y2=Node(Elem(ElemNo)%NodeNo(2))%Coord(2)
Elem(ElemNo)%L=sqrt((x1-x2)**2+(y1-y2)**2) ! 得到单元长度
cosA=(x2-x1)/Elem(ElemNo)%L
sinA=(y2-y1)/Elem(ElemNo)%L
Elem(ElemNo)%T=0.d0
Elem(ElemNo)%T(1,1)= cosA; Elem(ElemNo)%T(1,2)=sinA; ! 得到单元坐标转换矩阵
Elem(ElemNo)%T(2,1)=-sinA; Elem(ElemNo)%T(2,2)=cosA;
Elem(ElemNo)%T(3:4,3:4)=Elem(ElemNo)%T(1:2,1:2)
Elem(ElemNo)%B(1,1)=-1.d0/Elem(ElemNo)%L ! 得到单元几何矩阵
Elem(ElemNo)%B(1,2)=+1.d0/Elem(ElemNo)%L
Elem(ElemNo)%EPS=0.d0; Elem(ElemNo)%SIG=0.d0;
Elem(ElemNo)%dEPS=0.d0; Elem(ElemNo)%STRESS=0.d0;
Elem(ElemNo)%fy_temp=Elem(ElemNo)%fy; Elem(ElemNo)%Et=Elem(ElemNo)%E
Elem(ElemNo)%R=0.d0;
end do
return
end subroutine Intial_Elem
end module TypDef

module Main_Iteration ! 主迭代模块
use TypDef
implicit none

contains

subroutine Iteration(GValue, Node, Elem, Load, Support) !主迭代程序
type(typ_GValue) :: GValue
type(typ_Node) :: Node(:)
type(typ_Truss) :: Elem(:)
type(typ_Load) :: Load(:)
type(typ_Support) :: Support(:)
real(rkind) :: GLoad(2*GValue%NNode), GDisp(2*GValue%NNode) ! 总荷载向量，总位移向量
real(rkind) :: GK(2*GValue%NNode,2*GValue%NNode) ! 总刚度矩阵
real(rkind) :: RLoad(2*GValue%NNode) ! 节点反力
integer (ikind) :: NCycle ! 迭代次数
integer (ikind) :: i
NCycle=1
do
print *, " 第",NCycle,"次迭代"

call Get_GLoad(GValue,Load, GLoad) ! 得到当前荷载
call Get_RLoad(GValue, Elem, RLoad) ! 得到节点反力
GLoad=GLoad-RLoad ! 得到节点不平衡力
do i=1, GValue%NElem; call Get_EK(Elem(I)); end do ! 得到单元刚度矩阵
call Get_GK (GValue, Elem, GK) ! 得到整体刚度矩阵
call Get_Support(GValue, GK, GLoad, Support) ! 组装支座信息
call Check_Error(GValue, Load, GLoad) ! 检查是否收敛
print *, " 当前误差为",GValue%Error, "误差容限为", GValue%Tol
if(GValue%Error<=GValue%Tol) then
print *, " 当前步迭代收敛"
Elem(:)%EPS=Elem(:)%EPS+Elem(:)%dEPS ! 更新单元应变
Elem(:)%SIG=Elem(:)%STRESS ! 更新单元应力
Elem(:)%fy=Elem(:)%fy_temp ! 更新屈服应力
do i=1,size(Node)
Node(i)%Disp=Node(i)%Disp+Node(i)%dDisp ! 更新节点位移
Node(i)%dDisp=0.d0 ! 节点位移增量清零
end do
exit ! 退出迭代
end if
GK=matinv(GK) ! 刚度矩阵求逆
GDisp=matmul(GK,GLoad) ! 得到位移增量
call Get_dDisp(GValue, Node, GDisp) ! 得到位移增量
do i=1, GValue%NElem; call Get_EPS(Node,Elem(I)); end do ! 得到应变增量
do i=1, GValue%NElem; call Get_Stress(Elem(I)); end do ! 得到当前应力
do i=1, GValue%NElem; call Get_R(Elem(I)); end do ! 得到节点反力
NCycle=NCycle+1
end do
return
end subroutine Iteration

subroutine Check_Error(GValue, Load, GLoad) ! 检查是否收敛，采用无穷范数 e<|R|/|Load|
type(typ_GValue) :: GValue
type(typ_Load) :: Load(:)
real(rkind) :: GLoad(:)
real(rkind) :: x,y
integer (ikind) :: i
x=maxval(abs(GLoad)) ! 得到不平衡力的最大值
y=0.d0
do i=1,size(Load)
y=max(y,maxval(abs(Load(i)%Value))/real(GValue%NStep))
end do
GValue%Error=x/y ! 得到相对误差
return
end subroutine

subroutine Get_Support(GValue, GK, GLoad, Support) ! 采用化零置一法集成支座信息
type(typ_GValue) :: GValue
type(typ_Support) :: Support(:)
real(rkind) :: GK(:,:),GLoad(:)
integer(ikind) :: i,j
do i=1,size(Support)
j=Support(i)%NodeNo*2-2+Support(i)%DOF ! 得到支座对应的自由度编号
GK(j,:)=0.d0
GK(:,j)=0.d0
GK(j,j)=1;
GLoad(j)=0.d0
end do
return
end subroutine Get_Support

subroutine Get_GK (GValue, Elem, GK) ! 组装整体刚度矩阵
type(typ_GValue) :: GValue
type(typ_Truss) :: Elem(:)
real(rkind) :: GK(:,:)
integer(ikind) :: i,j,k,N(4)
GK=0.d0;
do i=1, size(Elem)
N=(/Elem(i)%NodeNo(1)*2-1, Elem(i)%NodeNo(1)*2, Elem(i)%NodeNo(2)*2-1, &
Elem(i)%NodeNo(2)*2/)
do j=1,4
do k=1,4
GK(N(j),N(k))=GK(N(j),N(k))+Elem(i)%K(j,k)
end do
end do
end do
return
end subroutine Get_GK

subroutine Get_RLoad(GValue, Elem, RLoad) ! 得到结构反力
type(typ_GValue) :: GValue
type(typ_Truss) :: Elem(:)
real(rkind) :: RLoad(:)
real(rkind) :: R(4)
integer (ikind) :: i
RLoad=0.d0
do i=1,size(Elem)
R=Elem(i)%R
RLoad(Elem(i)%NodeNo(1)*2-1)=RLoad(Elem(i)%NodeNo(1)*2-1)+R(1)
RLoad(Elem(i)%NodeNo(1)*2 )=RLoad(Elem(i)%NodeNo(1)*2 )+R(2)
RLoad(Elem(i)%NodeNo(2)*2-1)=RLoad(Elem(i)%NodeNo(2)*2-1)+R(3)
RLoad(Elem(i)%NodeNo(2)*2 )=RLoad(Elem(i)%NodeNo(2)*2 )+R(4)
end do
return
end subroutine Get_RLoad

subroutine Get_GLoad(GValue,Load, GLoad) ! 得到总体外荷载
type(typ_GValue) :: GValue
type(typ_Load) :: Load(:)
real(rkind) :: GLoad(:)
integer (ikind) :: i
GLoad=0.d0
do i=1, size(Load)
GLoad(Load(i)%NodeNo*2-1)=GLoad(Load(i)%NodeNo*2-1)+Load(i)%Value(1)
GLoad(Load(i)%NodeNo*2 )=GLoad(Load(i)%NodeNo*2 )+Load(i)%Value(2)
end do
GLoad=GLoad*real(GValue%CurrentStep)/real(GValue%NStep)
return
end subroutine Get_GLoad

function matinv(A) result (B) ! 对矩阵求逆
real(rkind) ,intent (in)::A(:,:)
real(rkind) , pointer ::B(:,:)
integer(ikind):: N,I,J,K; real(rkind)::D,T
integer(ikind), allocatable::IS(:),JS(:)
N=size(A,dim=2); allocate(B(N,N))
allocate(IS(N));allocate(JS(N))
B=A
do K=1,N
D=0.0D0
do I=K,N
do J=K,N
if(abs(B(I,J))>D) then
D=abs(B(I,J)); IS(K)=I; JS(K)=J
end if
end do
end do
do J=1,N
T=B(K,J); B(K,J)=B(IS(K),J); B(IS(K),J)=T
end do
do I=1,N
T=B(I,K); B(I,K)=B(I,JS(K)); B(I,JS(K))=T
end do
B(K,K)=1.d0/B(K,K);
do J=1,N; if(J.NE.K) B(K,J)=B(K,J)*B(K,K); end do
do I=1,N
if(I.NE.K) then
do J=1,N; if(J.NE.K) B(I,J)=B(I,J)-B(I,K)*B(K,J); end do
end if
end do
do I=1,N; if(I.NE.K) B(I,K)=-B(I,K)*B(K,K); end do
end do
do K=N,1,-1
do J=1,N
T=B(K,J); B(K,J)=B(JS(K),J); B(JS(K),J)=T
end do
do I=1,N
T=B(I,K); B(I,K)=B(I,IS(K)); B(I,IS(K))=T
end do
end do
return
end function matinv

end module Main_Iteration

module Data_Input
use TypDef
implicit none

contains

subroutine DataInput(GValue, Node, Elem, Load, Support) ! 读入数据
type(typ_GValue) :: GValue
type(typ_Node),pointer :: Node(:)
type(typ_Truss),pointer :: Elem(:)
type(typ_Load),pointer :: Load(:)
type(typ_Support),pointer :: Support(:)
integer (ikind) :: i,j
open(55,file='input.txt')
read(55,*) GValue%NStep, GValue%Tol ! 读入荷载步数，误差容限
read(55,*) GValue%NNode ! 读入节点数目
allocate(Node(GValue%NNode))
do i=1, GValue%NNode
read(55,*) j, Node(i)%Coord ! 读入节点坐标
Node(i)%Disp=0.d0; Node(i)%dDisp=0.d0
end do

read(55,*) GValue%NElem ! 读入单元数目
allocate(Elem(GValue%NElem))
do i=1, GValue%NElem
read(55,*) j, Elem(i)%NodeNo(1), Elem(i)%NodeNo(2), Elem(i)%A, &
Elem(i)%E, Elem(i)%fy, Elem(i)%Eh
! 读入单元节点编号，单元截面积，初始弹性模量，屈服强度，硬化模量。
end do
call Intial_Elem(Node,Elem) !赋予单元初始条件

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

read(55,*) GValue%NSupport ! 读入支座数
allocate(Support(GValue%NSupport))
do i=1, GValue%NSupport
read(55,*) j, Support(i)%NodeNo, Support(i)%DOF ! 读入节点，支座约束自由度
end do

close(55)
return
end subroutine DataInput

end module Data_Input

module Data_Output
use TypDef
implicit none

contains
subroutine DataOutput(GValue, Node, Elem, Load, Support) ! 计算结果输出模块
type(typ_GValue) :: GValue
type(typ_Node),pointer :: Node(:)
type(typ_Truss),pointer :: Elem(:)
type(typ_Load),pointer :: Load(:)
type(typ_Support),pointer :: Support(:)
integer (ikind) :: i
open(55,file="dispout.txt")
do i=1, size(Node)
write(55,'(I3,2G13.5)') i, Node(i)%Disp(1:2) ! 输出最终节点位移
end do
close(55)
open(55,file="Elemout.txt")
do i=1,size(Elem)
write(55,'(I3,2G13.5)') i, Elem(i)%EPS, Elem(i)%SIG ! 输出最终单元应变，应力
end do
close(55)
return
end subroutine DataOutput
end module Data_Output

program Main
use TypDef
use Main_Iteration
use Data_Input
use Data_Output
type(typ_GValue) :: GValue
type(typ_Node),pointer :: Node(:)
type(typ_Truss),pointer :: Elem(:)
type(typ_Load),pointer :: Load(:)
type(typ_Support),pointer :: Support(:)
integer (ikind) :: Inc ! 增量步

call DataInput(GValue, Node, Elem, Load, Support) ! 读入数据
do Inc=1, GValue%NStep
print *,"****************************************************"
print *,"当前增量步", inc
GValue%CurrentStep=Inc
call Iteration(GValue, Node, Elem, Load, Support) ! 进行迭代
end do
call DataOutput(GValue, Node, Elem, Load, Support) ! 输出结果
stop
end program