! ********************************************************************
! 一个牛顿－拉弗逊法求解非线性问题迭代过程的简单程序例子
! Example source code for Newton-Raphson Method in Non-linear Problems
! 可以进行荷载步长自动调整
! Load step can be automatically adjusted
! 读者可参考该模块加入相应的子程序编制自己的非线性计算程序
! 作 者：陆新征
! Author: Lu Xinzheng
! 指导教师：江见鲸
! Supervisor: Jiang Jianjing
! 清华大学土木工程系
! Dept. Civil Engrg. of Tsinghua University
! last revised: 2002.12.
!*********************************************************************

module Main_Iteration

use Lxz_Tools ! 通用子程序模块
use TypeDef ! 数据结构定义模块
use Data_Input ! 数据输入模块
use Data_Output ! 数据输出模块
use Elem_Prop ! 单元属性模块
use MatSolve ! 矩阵求解模块
implicit none

contains

subroutine RC2D_Main()
type(typ_GValue) :: GValue !总体控制变量数据结构数组
type(typ_Node),pointer :: Node(:) !节点数据结构数组
type(typ_Elem),pointer :: Elem(:) !混凝土单元数据结构数组
type(typ_Rebar),pointer :: Rebar(:) !钢筋单元数据结构数组
type(typ_Load),pointer :: Load(:) !荷载数据结构数组
type(typ_Load),pointer :: Initial_Load(:) !初始荷载
type(typ_Material),pointer :: Material(:) !材料数据结构数组
type(typ_Support),pointer :: Support(:) !支座数据结构数组

call DataInput(GValue,Node,Elem,Rebar,Material,Load,Support,&
Initial_Load) !读入数据

call Iteration(GValue,Node,Elem,Rebar,Load,Support,Initial_Load) !迭代核心程序

call DataOutput(GValue,Node,Elem,Rebar,Material,Load,Support) !输出数据

write(*,*) "计算成功结束"
read(*,*)
deallocate(node)
deallocate(Elem)
deallocate(Rebar)
deallocate(Load)
deallocate(Support)
return
end subroutine

subroutine Iteration(GValue0,Node0,Elem0,Rebar0,Load0,Support0,Initial_Load0)
type(typ_GValue) :: GValue0
type(typ_Node) :: Node0(:)
type(typ_Elem) :: Elem0(:)
type(typ_Rebar) :: Rebar0(:)
type(typ_Load):: Load0(:)
type(typ_Load):: Initial_Load0(:)
type(typ_Support) :: Support0(:)
type(typ_GValue) :: GValue
type(typ_Node) :: Node(GValue0%NNode)
type(typ_Elem) :: Elem(GValue0%NElem)
type(typ_Rebar) :: Rebar(GValue0%NRebar)
type(typ_Load):: Load(GValue0%NLoad)
type(typ_Load):: Initial_Load(GValue0%NInitial_Load)
type(typ_Support) :: Support(GValue0%NSupport)
Integer(ikind) :: I,J,II,III

GValue=GValue0
GValue%FinishedStep=0

II=1 !开始计算加载步数

LOOP_A: do
write(*,*) "======================================================="
write(*,*) 'Step',II,"总加载步数",GValue%NStep
write(*,*) '完成',(GValue%FinishedStep)/real(GValue%NStep)*100.0,'%'

!对数据进行备份，在荷载步调整操作时恢复上一步的参数
do I=1,GValue%NNode; Node(I)=Node0(I); end do
do I=1,GValue%NElem; Elem(I)=Elem0(I); end do
do I=1,GValue%NRebar; Rebar(I)=Rebar0(I); end do
do I=1,Gvalue%NLoad; Load(I)=Load0(I); end do
do I=1,Gvalue%NInitial_Load; Initial_Load(I)=Initial_Load0(I); end do
do I=1,GValue%NSupport; Support(I)=Support0(I); end do

GValue%GError=1.0 !初始误差设定为1.0

!=============================================================
GValue%ForError=0.0d0 !前一次迭代误差清零

III=1
LOOP_B: do
write(*,*) 'Iteration Number:', III
call Get_Elem_B(GValue,Elem,Node) !计算形函数
call Get_Elem_D(GValue,Elem) !计算混凝土单元材料本构矩阵
call Get_Elem_K(GValue,Elem) !计算混凝土单元刚度矩阵
call Get_Rebar_K(GValue,Rebar,Node) !计算钢筋单元刚度矩阵
call Get_Elem_F(GValue,Elem) !计算混凝土单元荷载向量
call Get_Rebar_F(GValue,Rebar) !计算钢筋单元荷载向量
!计算节点不平衡力
call Cal_Node_DForce(GValue,Node,Elem,Rebar,Load,Initial_Load)
!计算荷载误差
call Get_Error_F(GValue, Node, Load, Initial_Load,support)

if(III.ne.1) then !开始迭代
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! call Get_Error_F(GValue, Node)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
write(*,*) "当前迭代误差:",GValue%GError
if(GValue%GError<=GValue%ErrorT) then !误差小于误差容限，进行下一步计算
!!!!!!!!!!!!!!!!!!!!!
do I=1,GValue%NElem
!更新相应的参数
do J=1,4
Elem(I)%Concrete(J)%EPS=Elem(I)%Concrete(J)%STRAIN
Elem(I)%Concrete(J)%SIG=Elem(I)%Concrete(J)%Stress
Elem(I)%SRebar(J)%EPS=Elem(I)%SRebar(J)%STRAIN
Elem(I)%SRebar(J)%SIG=Elem(I)%SRebar(J)%Stress
end do
end do
do I=1, GValue%NRebar
Rebar(I)%SIG=Rebar(I)%Stress
Rebar(I)%EPS=Rebar(I)%Strain
end do
!把新的参数赋给保留参数
do I=1,GValue%NNode; Node0(I)=Node(I); end do
do I=1,GValue%NElem; Elem0(I)=Elem(I); end do
do I=1,GValue%NRebar; Rebar0(I)=Rebar(I); end do
do I=1,Gvalue%NLoad; Load0(I)=Load(I); end do
do I=1,GValue%NSupport; Support0(I)=Support(I); end do
write(*,*) "计算收敛,进行下一步计算"
II=II+1
GValue%FinishedStep=GValue%FinishedStep+1 !完成一步计算
exit LOOP_B
end if
end if
!得到前3此误差修正
if(III==2) GValue%ForError(1)=GValue%GError
if(III==3) GValue%ForError(2)=GValue%GError
if(III==4) GValue%ForError(3)=GValue%GError
if(III>4) then
GValue%ForError(1)=GValue%ForError(2)
GValue%ForError(2)=GValue%ForError(3)
GValue%ForError(3)=GValue%GError
end if
if(III>=30 .and. GValue%ForError(1)<GValue%ForError(2) .and. &
GValue%ForError(2)<GValue%ForError(3).or.&
GValue%GError>1.0d6) then !三次误差持续增大
!荷载折半
GValue%NStep=GValue%NStep*2

if(GValue%NStep>GValue%MaxStep) then
write(*,*) "最大分割步数,计算不收敛"
stop
end if

GValue%FinishedStep=GValue%FinishedStep*2
write(*,*) "计算收敛困难,步长减半"
exit LOOP_B
end if
!得到整体刚度矩阵并求解
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
call Get_GK(GValue,Node,Elem,Rebar,Load,Support,Initial_Load)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
III=III+1
end do LOOP_B !for III
!计算成功结束
if(GValue%FinishedStep>=GValue%NStep) then
exit LOOP_A
end if
end do LOOP_A !for II
return
end subroutine Iteration

end module