01变带宽总刚集成及求解简单程序例子.docVIP

01 变带宽总刚集成及求解简单程序例子 C ============================================ C 变带宽总刚集成及求解简单程序例子 C LU分解法，in-core only C 主要参考文献： C 1. 袁驷 程序结构力学 C 2. 江见鲸 陆新征 叶列平 混凝土结构有限元分析 C ============================================ 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 ! 罚函数法确定边