module TypeDef !定義混凝土材性模塊
type :: typ_Concrete
real*8 fc, ft,E0, ENU, EPS_Crush !抗壓強度+，抗拉強度+，初始彈性模量，初始泊松比，壓碎應變-
real*8 Es, ENUs !割線模量，割線泊松比
real*8 T(6,6) !坐標轉換矩陣
integer NCrack (3), Pre_NCrack(3), Pre_inc, Pre_incsub; !開裂記錄，前次迭代開裂記錄，前次增量步，前次增量子步
real*8 SIG(6), EPS(6),dEPS(6); !開始時應力，應變，應變增量
real*8 StressP(6), StrainP(6); !主應力，主應變
real*8 Stress(6), Strain(6) !結束時應力，應變
real*8 Beta,Pre_Beta !非線性指標，前次迭代非線性指標
real*8 D(6,6), Dela(6,6), Ds(6,6) !剛度矩陣，彈性剛度矩陣，割線剛度矩陣
end type typ_Concrete
end module TypeDef
module My_MOD !開辟公共變量空間
use TypeDef
type(typ_Concrete) :: My_Con(1000,8) !定義混凝土數組
end module


subroutine Get_DS(D,G,E,DE,S,TEMP0,
1 DTEMP,NGENS,N,NN,KC,MATS,NDI,NSHEAR,inc,incsub,ncycle)
! D(6x6) 迭代本構矩陣(out)
! G(6) 由于狀態改變引起的應力變化，不用(out)
! E(6) 開始時刻的應變(in)
! DE(6) 應變增量(in)
! S(6) 開始時刻的應變(in & out)
! Temp0 溫度(in)
! DTEMP 溫度變化(in)
! NGENS 應變維數(in)
! N(2) 單元編號(in)
! NN 積分點編號(in)
! KC 層號(in)
! MATS 材料編號(in)
! NDI 正應力維數(in)
! NSHEAR 剪應力維數(in)
! inc 當前增量步(in)
! incsub 當前增量子步(in)
! ncycle 當前循環數(in)

use IMSL !引用IMSL函數庫
use typedef
use My_Mod
implicit none
integer :: ngens,nn,kc,mats,ndi,nshear,inc,incsub,ncycle
real*8 :: e(ngens),de(ngens),temp0(1),dtemp(1),g(ngens)
1 ,d(ngens,ngens),s(ngens)

integer :: n(2)

type(typ_concrete) :: C
real*8 Beta1,strain_m
real*8 s_m,J2,J3,r,sita,TempA, TempB, TempC;
integer NSubStep !子步積分步數
integer I,J,K1,K2


C=My_Con(n(1),nn) !得到內存中保留的數據

c ----------------------------------------------------
c 初次計算，清零并賦值
if(inc==0.and.incsub==0.and.ncycle==0) then
open(77,file='debug.txt')
write(77,*)
close(77)

C%fc=30.; C%ft=3.; C%E0=30e3; C%ENU=0.18;
C%EPS_Crush=-0.0033;
C%T=0.; C%NCrack=0; C%Pre_NCrack=0;
C%Beta=0; C%Pre_Beta=0;
C%Pre_inc=0;
C%Pre_incsub=0;
C%SIG=0.; C%EPS=0.; C%Stress=0.; C%Strain=0.;
end if
c ----------------------------------------------------
c 如果新的增量步開始，則更新相應變量
if(inc>C%Pre_inc .or. incsub>C%Pre_incsub) then
C%Pre_inc=inc; C%Pre_incsub=incsub
C%NCrack=C%Pre_NCrack; !修正裂縫狀態
C%Beta=C%Pre_Beta; !修正非線性指標狀態
! 判斷是否壓壞
strain_m=(C%EPS(1)+C%EPS(2)+C%EPS(3))/3.
if(Strain_m>0.) Strain_m=0.
if(minval(C%EPS(1:3)-Strain_m)<C%EPS_Crush) then
C%NCrack=100 !徹底破壞
end if
end if

c 數據賦值
open(77,file='debug.txt',position='append')

C%SIG=s; C%EPS=e; C%dEPS=de;
C%Pre_NCrack=C%NCrack
C%Pre_Beta=C%Beta
NSubStep=4
c ------------------------------------
c 計算彈性矩陣
C%Dela=0.
do K1=1, 3
do K2=1, 3
C%Dela(K2,K1)=C%ENU
end do
C%Dela(K1,K1)=1.-C%ENU
END do
do K1=4,6
C%Dela(K1,K1)=(1.-2.*C%ENU)*0.5
end do
C%Dela=C%Dela*C%E0/(1.+C%ENU)/(1.-2.*C%ENU)
c --------------------------------------
c 如果已經壓碎，應力清零，剛度為很小值，結束計算
if(maxval(C%NCrack)==100) then
C%D=0.0001*C%Dela
C%SIG=0.;
C%Stress=0.;
s=0.
return
end if

C 計算主應力和割線剛度
C%Stress=C%SIG
do I=1, NSubStep

s_m=(C%Stress(1)+C%Stress(2)+C%Stress(3))/3. !計算平均應力
s(1:3)=C%Stress(1:3)-s_m
s(4:6)=C%Stress(4:6) !計算應力偏量
J2=-s(1)*s(2)-s(2)*s(3)-s(3)*s(1)+s(4)**2+s(5)*2+s(6)**2 !計算J2
J3=s(1)*s(2)*s(3)+2.*s(4)*s(5)*s(6) ! 計算J3
1 -s(1)*s(5)**2-s(2)*s(6)**2-s(3)*s(4)**2
r=sqrt(4.*J2/3.)
if(r.ne.0.) then
sita=acos(4.*J3/r**3)/3.
else
sita=0.
end if

if(maxval(abs(C%Pre_NCrack))==0) then !沒有裂縫
call Get_T_Matrix(C%SIG,C%T) !計算坐標轉換矩陣
end if

C%StressP=matmul(transpose(C%T),C%Stress); !計算主應力
C%StrainP=matmul(transpose(C%T),C%Strain); !計算主應變

if(C%StressP(1)<0.05*C%fc.and.
1 maxval(abs(C%Pre_NCrack))==0) then !沒有裂縫
TempA=1.2856/C%fc**2;
TempB=(1.4268+10.2551*cos(sita))/C%fc;
TempC=3.2128*s_m*3./C%fc-1.;
Beta1=-TempB+sqrt(TempB**2-4.*TempA*TempC)
Beta1=Beta1/2./TempA
Beta1=sqrt(J2)/Beta1 !根據江見鯨模型求解Beta

if(Beta1>C%Beta) then !Beta應該始終增大(對于全量模型)
C%Pre_Beta=Beta1
else
Beta1=C%Beta
end if
! 計算割線剛度和泊松比
if(Beta1>1.) Beta1=1.
C%Es=C%E0/2.*(1.+sqrt(1.-Beta1))
if(Beta1<0.8) C%ENUs=C%ENU
if(Beta1.ge.0.8) C%ENUs=0.42-(0.42-C%ENU)*
1 sqrt(1.-((Beta1-.8)/.2)**2)
!計算割線剛度矩陣
C%Ds=0.
do K1=1, 3
do K2=1, 3
C%Ds(K2,K1)=C%ENUs
end do
C%Ds(K1,K1)=1.-C%ENUs
end do
do K1=4,6
C%Ds(K1,K1)=(1.-2.*C%ENUs)*0.5
end do
C%Ds=C%Ds*C%Es/(1.+C%ENUs)/(1.-2.*C%ENUs)
else !如果處于開裂控制區
C%Ds=C%Dela
do K1=1,3
if(C%StressP(K1)>C%ft .OR. C%Pre_NCrack(K1)>0) then !按開裂處理
C%Pre_NCrack(K1)=1;
Call Crack_Open(C,K1) !計算開裂矩陣
end if
end do

C%Ds=matmul(C%T,matmul(C%Ds,transpose(C%T))); !計算割線剛度矩陣

end if
if(Beta1<0.99999d0) then !如果沒有達到極限應力
C%Strain=C%EPS+C%dEPS*real(I/NSubStep);
C%Stress=matmul(C%Ds,C%Strain);

else !達到極限應力后應力不變
C%Stress=C%SIG
C%Ds=1.d-6*C%Ds
end if
end do

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c 設置迭代剛度矩陣
D=c%Ds+1.d-6*C%Dela
s=C%Stress
My_Con(N(1),nn)=C;

c
c write(77,*) 'Step: ',inc, incsub, ncycle
c write(77,*) Beta1
c write(77,*) C%Pre_NCrack
c write(77,*) C%Stress(1:3)
c write(77,*)
close(77)
return
end subroutine

C 根據開裂修正剛度矩陣
subroutine Crack_Open(C,I0)
use typeDef
implicit none
type (typ_Concrete) :: C
integer,intent(in) :: I0
real*8 FACTOR,ECRA1,ECRA2,ECRA12

C%Ds(I0,:)=0.
C%Ds(:,I0)=0.

ECRA1=C%StrainP(I0);
ECRA2=0.;
ECRA12=maxval(abs(C%StrainP(4:6)))
c 計算裂面剪力傳遞系數
call USHRET0 (FACTOR,ECRA1,ECRA2,ECRA12)
C%Ds(4,4)=FACTOR*C%Ds(4,4)
C%Ds(5,5)=FACTOR*C%Ds(5,5)
C%Ds(6,6)=FACTOR*C%Ds(6,6)

return
end subroutine


c 計算裂面剪力傳遞系數
SUBROUTINE USHRET0 (FACTOR,ECRA1,ECRA2,ECRA12)
IMPLICIT REAL *8 (A-H, O-Z)
factor=0.4;
RETURN
END

c 計算主應力及坐標轉換矩陣
subroutine Get_T_Matrix(olds,T)
use IMSL
implicit none
real*8 olds(6), T(6,6)
real*8 SIG(3,3),eval_r(3), EVEC(3,3)
real*8 l(3),m(3),n(3)
real*8 SIGP(3)
integer I
SIG(1,1)=olds(1)
SIG(2,2)=olds(2)
SIG(3,3)=olds(3)
SIG(1,2)=olds(4); SIG(2,1)=olds(4);
SIG(2,3)=olds(5); SIG(3,2)=olds(5);
SIG(1,3)=olds(6); SIG(3,1)=olds(6);
call DEVCSF(3, SIG, 3, EVAL, EVEC, 3)
SIGP(1)=maxval(EVAL)
do I=1,3
if(eval_r(I)==SIGP(1)) then
l=EVec(:,I)
eval_r(I)=minval(EVAL)-10;
exit
end if
end do
SIGP(2)=maxval(EVAL)
do I=1,3
if(eval_r(I)==SIGP(2)) then
m=EVec(:,I)
eval_r(I)=minval(EVAL)-10;
exit
end if
end do
SIGP(3)=maxval(EVAL)
do I=1,3
if(eval_r(I)==SIGP(3)) then
n=EVec(:,I)
eval_r(I)=minval(EVAL)-10;
exit
end if
end do
do I=1,3
T(I,:)=(/l(I)**2,m(I)**2,n(I)**2,l(I)*m(I),
1 m(I)*n(I),n(I)*l(I)/)
end do
T(4,:)=(/2.d0*l(1)*l(2),2.d0*m(1)*m(2),2.d0*n(1)*n(2),
1 l(1)*m(2)+l(2)*m(1),m(1)*n(2)+m(2)*n(1),n(1)*l(2)+n(2)*l(1)/)
T(5,:)=(/2.d0*l(2)*l(3),2.d0*m(2)*m(3),2.d0*n(2)*n(3),
1 l(2)*m(3)+l(3)*m(2),m(2)*n(3)+m(3)*n(2),n(2)*l(3)+n(3)*l(2)/)
T(6,:)=(/2.d0*l(3)*l(1),2.d0*m(3)*m(1),2.d0*n(3)*n(1),
1 l(3)*m(1)+l(1)*m(3),m(3)*n(1)+m(1)*n(3),n(3)*l(1)+n(1)*l(3)/)
return
end subroutine


c marc 接口程序
SUBROUTINE HYPELA(D,G,E,DE,S,TEMP0,
1 DTEMP,NGENS,N,NN,KC,MATS,NDI,NSHEAR)
implicit real*8 (a-h,o-z)
INCLUDE '../common/concom' ! 通過concom模塊得到當前的計算步數
integer :: ngens,nn,kc,mats,ndi,nshear
real*8 :: e(1),de(1),temp0(*),dtemp(*),g(1),d(ngens,ngens),s(1)
integer :: n(2)
if(mats==1) then !如果材料編號是1
call Get_DS(D,G,E,DE,S,TEMP0,
1 DTEMP,NGENS,N,NN,KC,MATS,NDI,NSHEAR,inc,incsub,ncycle)
end if
return
end subroutine

c 后處理子程序
subroutine plotv(v,s,sp,etot,eplas,ecreep,t,m,nn,layer,ndi,
* nshear,jpltcd)
c* * * * * *
c
c select a variable contour plotting (user subroutine).
c
c v variable
c s (idss) stress array
c sp stresses in preferred direction
c etot total strain (generalized)
c eplas total plastic strain
c ecreep total creep strain
c t current temperature
c m(1) user element number
c m(2) internal element number
c nn integration point number
c layer layer number
c ndi (3) number of direct stress components
c nshear (3) number of shear stress components
c
c* * * * * *
use My_Mod
implicit real*8 (a-h,o-z) dp
dimension s(*),etot(*),eplas(*),ecreep(*),sp(*),m(2)
type(typ_Concrete) :: C
C=My_Con(m(1),nn)
c 后處理變量1：輸出是否開裂
if(jpltcd==1) then
if(C%NCrack(1).ne.0) then
v=1
else
v=0
end if
end if
c 后處理變量2-4：輸出裂縫狀態
do I=1,3
if(jpltcd==I+1) then
v=C%NCrack(I)
end if
end do
c 后處理變量5：輸出非線性指標
if(jpltcd==5) then
v=C%Beta
end if

c 后處理變量8：輸出裂縫數量
if(jpltcd==8) then
v=0.
do I=1,3
if(C%NCrack(I).ne.0) v=v+1
end do
end if

return
end