!============================================================================
Program RcConFea2005
!     Editor by Yiweiwen
!============================================================================
   use CoreCalculation
   implicit none
   integer(kind=ikind)::Nelem,Ntotv,Npoin,Nvfix,NSelem,NCelem,NFelem,Load_nodes,NLstep,istep
   integer(kind=skind)::Nnode,Ndofn,Nevab,Ngaus,type1,type2,type3,NSmat,NFmat,curvetype,LOCATECHAR
   integer(kind=skind),parameter::Ndime=2,Nstre=3

   logical if_file_exist,selfweight
   character(len=80):: filename,inputstr,refinestr
   real(kind=rkind) time_begin,time_end
   character(len=12)::Real_clock(3)
   integer(kind=ikind):: date_time(8)
   real(kind=rkind)::Fcu,Eo,Con_t,rzh,poisson !立方強度 初始彈性模量 單元厚度 容重 泊松比
   
!亂來變量   
   real(kind=rkind)::fc,Eiu,Etu,Ecu,Ec0,q 
   integer(kind=ikind)::i,j,k
   
   real(kind=rkind)::YL4(3),YL3(6),YL1(3)
!亂來變量   
    
   type(typ_Kcol),allocatable::Kcol(:)  !剛度
   type(Tptload),allocatable::ptload(:) !節點荷載

   integer(kind=ikind)::ESJtno(2),ESdof(4)  !鋼筋的節點好，自由度號
   integer(kind=ikind)::EFJtno(2),EFdof(4)   !粘結單元的
   real(kind=rkind)::ESStif(4,4),ESCoord(2,2),ESdis(4),ESF(4) !鋼筋的單元剛度，坐標，位移，力
   real(kind=rkind)::EFStif(4,4),EFCoord(2,2),EFdis(4),EFF(4) !粘結單元的   
   real(kind=rkind)::ECYB(3),ECYL(3)     !混凝土單元的應變和應力   
   real(kind=rkind)::T(3,3),Dmatx(3,3) !應力轉換矩陣和彈性矩陣

   real(kind=rkind),allocatable::G_Pcoord(:,:) !節點荷載矩陣
   real(kind=rkind),allocatable::G_load(:)     !荷載向量
   real(kind=rkind),allocatable::All_dis(:)    !總位移
   real(kind=rkind),allocatable::Ub_G_load(:)  !不平衡荷載向量
   integer(kind=ikind),allocatable::G_Pdof(:,:) !節點定位向量
  
   real(kind=rkind),allocatable::ECstif(:,:),ECcoord(:,:),ECdis(:) !混凝土的單元剛度，坐標，位移，力
   
   integer(kind=ikind),allocatable::Ecdof(:),ECJtno(:) !混凝土單元的節點好，自由度號

   integer(kind=ikind),allocatable::G_Ecdof(:,:),G_ECJtno(:,:) !混凝土單元整體的節點好，自由度號
   integer(kind=ikind),allocatable::G_ESdof(:,:),G_ESJtno(:,:)  !鋼筋的
   integer(kind=ikind),allocatable::G_EFdof(:,:),G_EFJtno(:,:) !粘結單元的

   integer(kind=skind),allocatable::SState(:),FState(:),CState(:,:) !鋼筋，粘結混凝土單元的狀態
   
   real(kind=rkind),allocatable::E(:,:,:),pu(:,:) !混凝土單元的E1，E2，和possion ratio
   real(kind=rkind),allocatable::diag(:),Bmatx(:,:),smatx(:,:) !一個是來分解方程，一個是應變矩陣

   integer(kind=skind),allocatable::SType(:),FType(:) !鋼筋和粘結各單元的材料組號向量
   real(kind=rkind),allocatable::steprop(:,:),feltprop(:,:) !鋼筋和粘結的材料向量

   real(kind=rkind),allocatable::SYL(:) !鋼筋的應力

   real(kind=rkind),allocatable::FYB(:,:),FYL(:),FForec(:,:),R_FYL(:)
   !粘結單元的應變，用前一步算出的應力，力，用當前步位移算出的應力
   real(kind=rkind),allocatable::CYB(:,:,:),C_main_YB(:,:,:),CYL(:,:,:),R_CYL(:,:,:)
   !混凝土單元的應變，主應變，用前一步算出的應力，用當前步位移算出的應力

   real(kind=rkind),allocatable::Fz(:,:,:),sic(:,:,:),siu(:,:,:)
   !存兩個方向的破壞應力，與破壞應力對于的應變，等效單向應變
 
 
!過程變量   
   real(kind=rkind)::area,thera,Es,Esl!混凝土單元面積，混凝土第一主應力的角度，鋼筋的彈性模量,粘結單元的滑移模量
   real(kind=rkind)::s,c,L !一般為角度的sin和cos,和桿的長度
   real(kind=rkind)::KH,KV !粘結單元的水平和垂直滑移模量
   real(kind=rkind)::tmp1,tmp2,tmp3,tmp4,tmp5  !亂來的變量
   real(kind=rkind)::s0,tmax  !粘結單元的H方向最大應力對應的位移，和最大應力
   integer(kind=skind)::Tx,Ty
   real(kind=rkind)::error,wuca  !容許誤差和計算過程中的誤差
   real(kind=rkind)::sl(3)  !真正的應力
   integer(kind=ikind)::numb !在每一個荷載步迭代的數目
   real(kind=rkind)::EFYB(2)  !在某一步中的粘結單元的應變增量
   real(kind=rkind)::angle
   integer(kind=ikind)::destroyNum,oldstate
   !integer(kind=ikind),allocatable::destroyNelem(:)
   !iteger(kind=ikind)::Dnelem

   
   

   
!program begin

   write(*,*)  
   write(*,"(a)") '!==========歡迎使用ReConFea2005程序,本程序由華南理工大學易偉文編寫==========!'
   write(*,"(a)") '!==========     RcConFea2005是一個鋼筋混凝土的非線性有限元程序   ===========!'
   write(*,*)
   write(*,"(a)") '請輸入數據文件的路徑：'
   filename="ConFea.txt"
   write(*,"(a)") 'ConFea.txt'
   read(*,"(a)") filename
   if (filename=='') filename="ConFea.txt"
   inquire(file=filename,Exist=if_file_exist)
   if (.not.if_file_exist) then
      write(*,*) "Error,file does not exist!"
	  stop
   end if
   open (unit=10 , file = filename , status = 'old' ,    action ='read')
   write(*,"(a)") '請輸入輸出文件的路徑：'
   filename="ConFeaout.txt"
   write(*,"(a)") 'ConFeaout.txt'
   read(*,"(a)") filename
   if (filename=='') filename="ConFeaout.txt"   
   open (unit=11 , file = filename , status = 'replace', action='readwrite') 
   
   open (unit=111 , file = "state.txt",status = 'replace', action='readwrite') 

   call date_and_time(real_clock(1),real_clock(2),real_clock(3),Date_time)
   call cpu_time(time_begin)
   
   READ(10,"(a)") inputstr
   LOCATECHAR=INDEX(inputstr,'=')
   LOCATECHAR=LOCATECHAR+1
   refinestr=inputstr(LOCATECHAR:80)
   READ(refinestr,*) filename

   write(11,"(a)") "============================================================================================"
   write(11,"(a)") "這個是ReConFea2005的結果輸出文件"
   write(11,"(a)") "============================================================================================"   
   write(11,"(A,14x,i5,A,i2.2,a,i2.2)") "Today is:",DAte_time(1),'-',date_time(2),'-',date_time(3)
   write(11,"(3(a,i2),a,i3.3)") 'The beginning time is:  ',Date_time(5),':',date_time(6),':',Date_time(7),'.',date_time(8)
   write(11,*)
   write(11,"(a)") "============================================================================================"  
   write(11,"(2a)") "本工程名稱：",filename
   write(11,"(a)") "============================================================================================"  

   READ(10,"(a)") inputstr
   LOCATECHAR=INDEX(inputstr,'=')
   LOCATECHAR=LOCATECHAR+1
   refinestr=inputstr(LOCATECHAR:80)
   READ(refinestr,*) type1
   if (type1==1) then
      write(11,"(a)") "下面進行的是彈性計算。"
   else
      write(11,"(a)") "下面進行的是彈塑性計算。"
   end if

   READ(10,"(a)") inputstr
   LOCATECHAR=INDEX(inputstr,'=')
   LOCATECHAR=LOCATECHAR+1
   refinestr=inputstr(LOCATECHAR:80)
   READ(refinestr,*) type2
   if (type2==1) then
      write(11,"(a)") "此工程為平面應力問題。"
   else
      write(11,"(a)") "此工程為平面應變問題。"
   end if

   READ(10,"(a)") inputstr
   LOCATECHAR=INDEX(inputstr,'=')
   LOCATECHAR=LOCATECHAR+1
   refinestr=inputstr(LOCATECHAR:80)
   READ(refinestr,*) type3
   selectcase(type3)
   case(1)
      write(11,"(a)") "本工程的混凝土單元使用T3單元，即常應變三角形單元。"
   case(2)
      write(11,"(a)") "本工程的混凝土單元使用Q4單元，即雙線性四邊形單元。"
   case(3)
      write(11,"(a)") "本工程的混凝土單元使用GC－Q6單元，即用廣義協調條件修正的Willon廣義協調四邊形單元。"
   case(4)
      write(11,"(a)") "本工程的混凝土單元使用GR12單元，即具有旋轉自由度的廣義協調四邊形單元。"
   case default
      write(11,"(a)") "本工程的混凝土單元使用GR12M單元，即在有旋轉自由度的基礎上引入泡沫位移的廣義協調四邊形單元。"
   end select  
   call confirm_parameter(type3,Nnode,Ngaus,Nevab,Ndofn)


   write(11,*)
   write(11,"(a)") "============================================================================================"  
   write(11,"(a)") "本工程總體信息："
   write(11,"(a)") "============================================================================================"  

   READ(10,"(a)") inputstr
   write(11,"(a)") inputstr
   LOCATECHAR=INDEX(inputstr,'=')
   LOCATECHAR=LOCATECHAR+1
   refinestr=inputstr(LOCATECHAR:80)
   READ(refinestr,*) NCelem

   READ(10,"(a)") inputstr
   write(11,"(a)") inputstr
   LOCATECHAR=INDEX(inputstr,'=')
   LOCATECHAR=LOCATECHAR+1
   refinestr=inputstr(LOCATECHAR:80)
   READ(refinestr,*) NSelem

   READ(10,"(a)") inputstr
   write(11,"(a)") inputstr
   LOCATECHAR=INDEX(inputstr,'=')
   LOCATECHAR=LOCATECHAR+1
   refinestr=inputstr(LOCATECHAR:80)
   READ(refinestr,*) NFelem
   Nelem=NCelem+NFelem+NSelem

   READ(10,"(a)") inputstr
   write(11,"(a)") inputstr
   LOCATECHAR=INDEX(inputstr,'=')
   LOCATECHAR=LOCATECHAR+1
   refinestr=inputstr(LOCATECHAR:80)
   READ(refinestr,*) Npoin

   READ(10,"(a)") inputstr
   write(11,"(a)") inputstr
   LOCATECHAR=INDEX(inputstr,'=')
   LOCATECHAR=LOCATECHAR+1
   refinestr=inputstr(LOCATECHAR:80)
   READ(refinestr,*) Nvfix

   READ(10,"(a)") inputstr
   write(11,"(a)") inputstr
   LOCATECHAR=INDEX(inputstr,'=')
   LOCATECHAR=LOCATECHAR+1
   refinestr=inputstr(LOCATECHAR:80)
   READ(refinestr,*) NLstep

   READ(10,"(a)") inputstr
   write(11,"(a)") inputstr
   LOCATECHAR=INDEX(inputstr,'=')
   LOCATECHAR=LOCATECHAR+1
   refinestr=inputstr(LOCATECHAR:80)
   READ(refinestr,*) selfweight

   READ(10,"(a)") inputstr
   write(11,"(a)") inputstr
   LOCATECHAR=INDEX(inputstr,'=')
   LOCATECHAR=LOCATECHAR+1
   refinestr=inputstr(LOCATECHAR:80)
   READ(refinestr,*) error

   READ(10,"(a)") inputstr
   write(11,"(a)") inputstr
   LOCATECHAR=INDEX(inputstr,'=')
   LOCATECHAR=LOCATECHAR+1
   refinestr=inputstr(LOCATECHAR:80)
   READ(refinestr,*) Curvetype

   READ(10,"(a)") inputstr
   write(11,"(a)") inputstr
   LOCATECHAR=INDEX(inputstr,'=')
   LOCATECHAR=LOCATECHAR+1
   refinestr=inputstr(LOCATECHAR:80)
   READ(refinestr,*) NSmat

   READ(10,"(a)") inputstr
   write(11,"(a)") inputstr
   LOCATECHAR=INDEX(inputstr,'=')
   LOCATECHAR=LOCATECHAR+1
   refinestr=inputstr(LOCATECHAR:80)
   READ(refinestr,*) NFmat

   READ(10,"(a)") inputstr
   write(11,"(a)") inputstr
   LOCATECHAR=INDEX(inputstr,'=')
   LOCATECHAR=LOCATECHAR+1
   refinestr=inputstr(LOCATECHAR:80)
   READ(refinestr,*) load_nodes

   allocate(G_Pcoord(2,Npoin),G_Pdof(Ndofn,Npoin))
   allocate(Ecstif(nevab,nevab),Eccoord(Ndime,Nnode),ECdis(nevab))
   allocate(Ecdof(nevab),ECJtno(Nnode))
   allocate(G_ESdof(4,NSelem),G_ESJtno(2,NSelem))
   allocate(G_EFdof(4,NFelem),G_EFJtno(2,NFelem))
   allocate(G_ECJtno(nnode,NCelem),G_Ecdof(nevab,NCelem))
   allocate(SState(NSelem),FState(NFelem))
   allocate(bmatx(3,nevab),smatx(3,nevab))
   allocate(SType(NSelem),FType(NFelem))
   allocate(steprop(3,NSmat),feltprop(4,NFmat))
   allocate(SYL(NSelem))
   allocate(FYB(2,NFelem),FYL(NFelem),R_FYL(NFelem),FForec(2,NFelem))

      
   if (type3==1) then       
	  allocate(CState(1,NCelem))	 
	  allocate(E(2,1,NCelem),pu(1,NCelem))
	  allocate(CYB(3,1,NCelem),C_main_YB(2,1,NCelem),CYL(6,1,NCelem),R_CYL(2,1,NCelem))
	  allocate(Fz(2,1,NCelem),sic(2,1,NCelem),siu(2,1,NCelem))
   else      
	  allocate(CState(4,NCelem))	  
	  allocate(E(2,4,NCelem),pu(4,NCelem))
	  allocate(CYB(3,4,NCelem),C_main_YB(2,4,NCelem),CYL(6,4,NCelem),R_CYL(2,4,NCelem))
	   allocate(Fz(2,4,NCelem),sic(2,4,NCelem),siu(2,4,NCelem))
  end if
         
   read(10,"(a)") inputstr
   read(10,"(a)") inputstr
   read(10,*)  Fcu,Eo,Con_t,rzh,poisson
   
   
   write(11,*)
   write(11,"(a)") "============================================================================================"  
   write(11,"(a)") "以下是本工程混凝土材料的特征參數：" 
   write(11,"(a)") "============================================================================================"  
   write(11,"(a)") "材料號      抗壓強度        初始彈性模量    單元厚度        容重            泊松比"
   write(11,"(i2,4x,5e16.4)") 1,Fcu,Eo,Con_t,rzh,poisson

   read(10,"(a)") inputstr
   read(10,"(a)") inputstr
   read(10,*) steprop

   write(11,*)
   write(11,"(a)") "============================================================================================"  
   write(11,"(a)") "以下是本工程鋼筋材料的特征參數：" 
   write(11,"(a)") "============================================================================================"  
   write(11,"(a)") "材料號      彈性模量        鋼筋面積        設計強度"   
   
   do i=1,NSmat
      write(11,"(i2,4x,3e16.4)") i,steprop(:,i)
   end do

   read(10,"(a)") inputstr
   read(10,"(a)") inputstr
   read(10,*) feltprop

   write(11,*)
   write(11,"(a)") "============================================================================================"  
   write(11,"(a)") "以下是本工程粘結單元的特征參數：" 
   write(11,"(a)") "============================================================================================"  
   write(11,"(a)") "材料號      鋼筋直徑和      控制長度        和x軸的角度     抗拉彈性模量"   
   
   do i=1,NFmat
      write(11,"(i2,4x,4e16.4))") i,feltprop(:,i)
   end do
   write(11,*)

   read(10,"(a)")inputstr
   write(11,"(a)") "============================================================================================"  
   write(11,"(a)") '鋼筋單元材料類型向量：'
   write(11,"(a)") "============================================================================================"  
   
   SType=1
   if(NSmat>1) then       
      read(10,*) SType	  
   end if
   do i=1,NSelem     
      write(11,"(i2)",advance='no') SType(i)      
   end do