! ********************************************************************
! 变带宽总刚集成及求解简单程序例子
! Source code example for in-core global matrix assemble and solve module
! LU分解法，in-core only
! 作 者：陆新征
! Author: Lu Xinzheng
! 指导教师：江见鲸
! Supervisor: Jiang Jianjing
! 清华大学土木工程系
! Dept. Civil Engrg. of Tsinghua University
! last revised: 2002.12.
!*********************************************************************
module MatSolve
use Lxz_Tools
use TypeDef
use Elem_Prop

implicit none

integer :: NGlbDOF !总自由度数

contains

subroutine Get_GK(GValue,Node,Elem,Rebar,Load,Support,Initial_Load)
type(typ_GValue) :: GValue ! 总体控制变量数据结构
type(typ_Node) :: Node(:) ! 节点性质数据结构数组
type(typ_Elem) :: Elem(:) ! 混凝土单元性质数据结构数组
type(typ_Rebar) :: Rebar(:) ! 钢筋单元性质数据结构数组
type(typ_Load) :: Load(:) ! 荷载数据结构数组
type(typ_Load) :: Initial_Load(:) ! 初始荷载数据结构数组
type(typ_Support) :: Support(:) ! 支座数据结构数组
type(typ_Kcol),allocatable :: Kcol(:) ! 变带宽总体刚度矩阵数组
real(rkind),allocatable :: GLoad(:), GDisp(:) ! 荷载，位移向量
integer(ikind) :: BandWidth(2*GValue%NNode) !带宽
real(rkind) :: Plenty ! 支座罚函数
integer(ikind) :: ELocVec(8)
integer(ikind) :: I,J,K
integer(ikind) :: MinDOFNum

NGlbDOF=2*GValue%NNode !计算总自由度数
Plenty=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,4
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,8
BandWidth(ELocVec(J))=min(MinDOFNum,BandWidth(ELocVec(J)))
end do
end do
do I=1,GValue%NRebar
do J=1,2
ELocVec(J*2-1)=2*Rebar(I)%NodeNo(J)-1
ELocVec(J*2 )=2*Rebar(I)%NodeNo(J)
end do
MinDOFNum=minval(ELocVec(1:4))
do J=1,4
BandWidth(ELocVec(J))=min(MinDOFNum,BandWidth(ELocVec(J)))
end do
end do

! 总刚矩阵分配内存
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,4
ELocVec(J*2-1)=2*Elem(I)%NodeNo(J)-1
ELocVec(J*2)=2*Elem(I)%NodeNo(J)
end do
do J=1,8
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
! 罚函数法确定边界条件
Plenty=max(Plenty,maxval(abs(Elem(I)%EK)))
end do

do I=1, GValue%NRebar
do J=1,2
ELocVec(J*2-1)=2*Rebar(I)%NodeNo(J)-1
ELocVec(J*2)=2*Rebar(I)%NodeNo(J)
end do
do J=1,4
do K=J,4
if(ELocVec(J)>ELocVec(K)) then
Kcol(ELocVec(J))%row(ELocVec(K))=&
Kcol(ELocVec(J))%row(ELocVec(K))+Rebar(I)%EK(J,K)
else
Kcol(ELocVec(K))%row(ELocVec(J))=&
Kcol(ELocVec(K))%row(ELocVec(J))+Rebar(I)%EK(J,K)
end if
end do
end do
Plenty=max(Plenty,maxval(abs(Rebar(I)%EK)))
end do
write(*,*) "完成总刚集成"

! 荷载向量集成
GLoad=0.d0

do I=1, GValue%NLoad
GLoad(Load(I)%NodeNo*2-1)=GLoad(Load(I)%NodeNo*2-1)&
+Load(I)%Value(1)/real(GValue%NStep)
GLoad(Load(I)%NodeNo*2 )=GLoad(Load(I)%NodeNo*2 )&
+Load(I)%Value(2)/real(GValue%NStep)
end do

write(*,*) "完成荷载集成"

! 支座向量集成
do I=1,GValue%NSupport
J=2*(Support(I)%NodeNo-1)+Support(I)%DOF
GLoad(J)=GLoad(J)+Support(I)%Value*Plenty*1.0D8
Kcol(J)%Row(J)=KCol(J)%Row(J)+Plenty*1.0d8
end do
write(*,*) "完成支座集成"

! 矩阵求解
call BandSolv(Kcol,Gload,GDisp)

! 得到节点位移增量
call Get_Node_dDisp(GValue,Node,GDisp)
! 得到单元应变增量
call Get_Elem_dEPS(GValue, Node, Elem)

! 关闭数组，释放内存
do I=1,NGlbDOF
deallocate(Kcol(I)%Row)
end do
deallocate(Kcol)
deallocate(GDisp)
deallocate(GLoad)
return
end subroutine Get_GK

! 得到节点位移增量
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)=GDisp(I*2-1)
Node(I)%dDisp(2)=GDisp(I*2 )
Node(I)%Disp=Node(I)%Disp+Node(I)%dDisp
end do
return
end subroutine Get_Node_dDisp

! 判断误差
subroutine Get_Error_F(GValue, Node,Load,Initial_Load,Support)
type(typ_GValue) :: GValue
type(typ_Node) :: Node(:)
type(typ_Load) :: Load(:)
type(typ_Load) :: Initial_Load(:)
type(typ_Support) :: Support(:)
real(rkind) :: External_Load(GValue%NNode*2)
integer :: I
real(rkind) :: dGForce(2*GValue%NNode)
real(rkind) :: Max_Force,Max_dForce
External_Load=0.
dGForce=0.

! 外力荷载向量集成
do I=1,GValue%NInitial_Load
External_Load(Initial_Load(I)%NodeNo*2-1)=&
External_Load(Initial_Load(I)%NodeNo*2-1)+Initial_Load(I)%Value(1)
External_Load(Initial_Load(I)%NodeNo*2 )=&
External_Load(Initial_Load(I)%NodeNo*2 )+Initial_Load(I)%Value(2)
end do
do I=1,GValue%NLoad
External_Load(Load(I)%NodeNo*2-1)=&
External_Load(Load(I)%NodeNo*2-1)+Load(I)%Value(1)*&
real(GValue%FinishedStep+1)/real(GValue%NStep)
External_Load(Load(I)%NodeNo*2 )=&
External_Load(Load(I)%NodeNo*2 )+Load(I)%Value(2)*&
real(GValue%FinishedStep+1)/real(GValue%NStep)
end do
! 节点不平衡力集成
do I=1, GValue%NNode
dGForce(I*2-1)=Node(I)%dForce(1)
dGForce(I*2 )=Node(I)%dForce(2)
end do
! 消去边界条件
do I=1, GValue%NSupport
dGForce((Support(I)%NodeNo-1)*2+Support%DOF)=0.
end do

! 力的零范数误差
if(GValue%ErrorNormal==0) then
Max_dForce=maxval(abs(dGForce))
Max_Force=maxval(abs(External_Load))
if(Max_Force.ne.0) then
GValue%GError=abs(Max_dForce/Max_Force)
else
GValue%GError=0.0
end if
end if

! 力的1范数误差
if(GValue%ErrorNormal==1) then
max_dForce=sum(abs(dGForce))
max_Force=sum(abs(External_Load))
if(Max_Force.ne.0) then
GValue%GError=abs(Max_dForce/Max_Force)
GValue%GError=GValue%GError/real(GValue%NNode)*&
real(GValue%NLoad)
else
GValue%GError=0.0
end if
end if

! 力的2范数误差
if(GValue%ErrorNormal==2) then
max_dForce=dot_product(dGForce,dGForce)
max_Force=dot_product(External_Load,External_Load)
if(Max_Force.ne.0) then
GValue%GError=abs(Max_dForce/Max_Force)
GValue%GError=GValue%GError/real(GValue%NNode**2)*&
real(GValue%NLoad**2)
else
GValue%GError=0.0
end if
end if
return
end subroutine Get_Error_F

! 矩阵求解

!---------------------------------------
subroutine BandSolv(Kcol,GLoad,Disp)
!---------------------------------------
type (typ_Kcol) :: Kcol(:);
real(rkind) :: GLoad(:),Disp(:);
integer :: row1,ncol,row,j,ie
real(rkind) :: diag(1:NGlbDOF),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
Disp(:)=GLoad(:)

return;
end subroutine BandSolv

end module