-
Notifications
You must be signed in to change notification settings - Fork 1
/
solvers.jl
218 lines (218 loc) · 10.1 KB
/
solvers.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
function solvers(u_inc::T, conf::T) where T<:Float64
## Initialization
m0(x) = (1.0.-x).^2 ./ (ka .- (ka .- 1.0).*(1.0 .- x).^2)
gd(x) = (1.0.-x).^2 ./ (ka .- ka.*(1.0 .- x).^2)
u = zeros(Float64,2*nnode)
u_old = zeros(Float64,2*nnode)
f_ext=zeros(Float64,2*nnode)
f_int=zeros(Float64,2*nnode)
f_int_old=zeros(Float64,2*nnode)
sigma=zeros(Float64,3,4*nel)
σc =zeros(Float64,3,4*nel)
CC = zeros(Float64,4*nel)
σd =zeros(Float64,3,4*nel)
σ_t = SharedArray{Float64,2}(12,nel) ## parallel computation of stress
σc_t = SharedArray{Float64,2}(12,nel)
σd_t = SharedArray{Float64,2}(3,4*nel)
epsilon = zeros(Float64,3,4*nel)
ψᵖ=zeros(Float64,4*nel)
εᵖ = zeros(Float64,3,4*nel)
εᵖ_old = zeros(Float64,3,4*nel)
δεᵖ = zeros(Float64,3,4*nel)
εᵖ_com = zeros(Float64,3,4*nel)
εᵖ_com_old = zeros(Float64,3,4*nel)
Hn1 = zeros(Float64,4,nel)
Hn1_old=copy(Hn1)
d1 = zeros(Float64,nnode).+d0; dg1=zeros(Float64,4*nel)
for iel = 1:nel
dg1[4*(iel-1)+1:4*iel]=Nd[:,4*(iel-1)+1:4*iel]*d1[dedofMat[iel,:]]
end
d1_old = deepcopy(d1)
dg1_old = deepcopy(dg1)
H0 = Initial_Hn(d1, dg1)
Hn1 = Hn1_comp!(H0,Hn1,ψᵖ,zeros(Float64,4*nel)) ## prescribed initial history field
err_storage = zeros(Float64,maxit,step_total+1) ##残差和收敛时间
time_storge = zeros(Float64,step_total+1) ##收敛时间
#Assemble whole Stiffness matrix
begin
@info "Formulating stiffness matrix takes"
@time KK=Kmatrix(element, Bu, detjacob,DK,iK,jK)
end
init_total = 1 ## step number in the first load inc
f_ext = FMat_ext(node,conf) ## external force
alldofs = 1:2size(node,1)
###
u_increment = u_inc
STOnit=zeros(Float64,step_total)
function monoli_initialStiff()
for stp=1:step_total
@info "step=: $stp"
r_ref = zeros(Float64,2*nnode)
u_ref = zeros(Float64,2*nnode)
du = zeros(Float64,2*nnode)
du[loaddofs] .= u_increment
## staggered iteration of u and d
err_d = 1.0; nit_d = 0
record = @timed while (err_d>1e-3) && (nit_d<maxit)
nit_d+=1
@info "nit_d=$nit_d"
# compute displacement field u
err=1; nit=0
δu = zeros(Float64,2*nnode)
rr = zeros(Float64,2*nnode)
while (err>tol) && (nit<maxit)
nit+=1
itd = 0
@info "-迭代步=: $nit"
if nit_d==1 && nit==1
rr[freedofs] .= KK[freedofs,loaddofs]*du[loaddofs] # .- (f_ext[freedofs]-f_int_old[freedofs])
r_ref .= rr
δu[freedofs] .= -KK[freedofs,freedofs]\(rr[freedofs])
u .= u_old .+ du .+ δu
u_ref .= du .+ δu
elseif nit>1
δu[freedofs] .= KK[freedofs,freedofs]\(rr[freedofs])
u .= u_old .+ δu
else
itd = 1
nothing
end
UU = u[edofMat]
# epsilon = zeros(Float64,3,4*nel)
for iel=1:nel
epsilon[:,4*(iel-1)+1:4*iel]=reshape(Bu[:,8*(iel-1)+1:8*iel]*UU[iel,:],3,4)
end
@sync @distributed for iel=1:nel
σ_t[:,iel] = blockdiag(sparse(reshape(DK[:,4*(iel-1)+1],3,3)),sparse(reshape(DK[:,4*(iel-1)+2],3,3)),
sparse(reshape(DK[:,4*(iel-1)+3],3,3)),sparse(reshape(DK[:,4*(iel-1)+4],3,3)))*
reshape((epsilon[:,4*(iel-1)+1:4*iel] .- εᵖ[:,4*(iel-1)+1:4*iel]),12,1)
σc_t[:,iel] = blockdiag(sparse(reshape(DK[:,4*(iel-1)+1],3,3)),sparse(reshape(DK[:,4*(iel-1)+2],3,3)),
sparse(reshape(DK[:,4*(iel-1)+3],3,3)),sparse(reshape(DK[:,4*(iel-1)+4],3,3)))*
reshape(εᵖ[:,4*(iel-1)+1:4*iel],12,1)
end
sigma .= reshape(sdata(σ_t),3,4*nel)
##local stress
σc .= sigma .- reshape(kron(gd.(dg1),ones(3)) .* sdata(σc_t[:]),3,4*nel)
###
CC .= sigma_P(σc,v,planetype) ## mean stress / microcrack open-closure criterion
idxTen = findall(CC.>=0.0)##拉伸断裂
idxCom = findall(CC.<0.0)##压剪断裂
if isempty(idxTen)==0 ##拉
# @info "拉伸"
εᵖ[:,idxTen] .= kron(1.0 .- m0.(dg1[idxTen]), ones(1,3))' .* (epsilon[:,idxTen])
end
# increment of plastic strian
δεᵖ .= plasticity!(σc,sigma,εᵖ,dg1,idxCom)
εᵖ[:,idxCom] .= εᵖ_old[:,idxCom] .+ δεᵖ[:,idxCom]
εᵖ_com .= εᵖ_com_old .+ δεᵖ # plastic strain due to friction
@sync @distributed for iel=1:nel
σ_t[:,iel] = blockdiag(sparse(reshape(DK[:,4*(iel-1)+1],3,3)),sparse(reshape(DK[:,4*(iel-1)+2],3,3)),
sparse(reshape(DK[:,4*(iel-1)+3],3,3)),sparse(reshape(DK[:,4*(iel-1)+4],3,3)))*
reshape(epsilon[:,4*(iel-1)+1:4*iel] .- εᵖ[:,4*(iel-1)+1:4*iel],12,1)
end
sigma .= reshape(sdata(σ_t),3,4*nel)
##
f_int = FMat(node,element,Bu,detjacob,sigma,edofMat)
# line search: secant line search
if abs((f_ext.-f_int)'*δu) > (0.8 * abs((f_ext.-f_int_old)' * δu)) && nit > 6
# line Search
# @info "-line Search"
nit_int = 0
LS = 0.0
LSS = 1.0
LSSS = 1.2
S = []
SS = []
SSS = []
id::Int64 = 0
err_ls=1
while err_ls > tol && nit_int < 15
nit_int += 1
LS = LSS
LSS = LSSS
#
SS = LSFunc(LSS,δu,u_old,εᵖ_old,d1_old,dg1_old,H0,Hn1_old,εᵖ_com_old,f_ext,r_ref)
S = LSFunc(LS,δu,u_old,εᵖ_old,d1_old,dg1_old,H0,Hn1_old,εᵖ_com_old,f_ext,r_ref)
LSSS = LSS - SS[1] / (SS[1] - S[1]) * (LSS - LS)
LSSS>30 || LSSS<-10 || abs(LSSS-LSS)<tol ? (@info "break" break) : nothing
LSSS>15.0 ? LSSS=15.0 : LSSS<0.1 ? LSSS=0.1 : nothing
##计算line search目标函数
# f_int .= SS[5]
err_ls = SS[1]/norm(r_ref)
@info "-line search err=$err_ls"
end
u .= SS[2]
sigma .= SS[3]
epsilon .= SS[4]
f_int .= SS[5]
d1 .= SS[6]
dg1 .= SS[7]
Hn1 .= SS[8]
εᵖ .= SS[9]
εᵖ_com .= SS[10]
CC .= SS[11]
end
rr[freedofs] .= f_ext[freedofs] .- f_int[freedofs]
# if itd == 1
# err = 1.0
# else
# err = abs((f_ext .- f_int)'*δu)/abs(r_ref' *u_ref)
# end
err = norm(rr[freedofs])/norm(r_ref)
@info "--err_u=$err"
f_int_old .= f_int
u_old .= u
d1_old .= d1
dg1_old .= dg1
Hn1_old .= Hn1
εᵖ_old .= εᵖ
εᵖ_com_old .= εᵖ_com
end
# compute damage field d
@sync @distributed for iel=1:nel
σc_t[:,iel] = blockdiag(sparse(reshape(DK[:,4*(iel-1)+1],3,3)),sparse(reshape(DK[:,4*(iel-1)+2],3,3)),
sparse(reshape(DK[:,4*(iel-1)+3],3,3)),sparse(reshape(DK[:,4*(iel-1)+4],3,3)))*
reshape(εᵖ[:,4*(iel-1)+1:4*iel],12,1)
end
σd .= reshape(sdata(σc_t),3,4*nel)
ψᵖ .= 1.0/2.0.*(σd[1,:].*εᵖ[1,:] .+ σd[2,:].*εᵖ[2,:] .+ σd[3,:].*εᵖ[3,:])
##计算d
Hn1 = Hn1_comp!(H0,Hn1,ψᵖ,CC)
d1, dg1 = integration_d(Hn1,d1,dg1)
# replace!(x->x>1.0 ? (1.0-k) : x, dg1)
# replace!(x->x<0.0 ? k : x, dg1)
err_d = norm(d1.-d1_old)/norm(d1)
# d1_old .= d1
# dg1_old .= dg1
@info "---err_d=$err_d"
err_storage[nit_d,stp] = err_d
end
# "Storing the output data takes:"
begin
Fload1[stp+1] = sum(f_int[loaddofs])/(size(loaddofs,1)-1)
Uload1[stp+1] = sum(u[loaddofs])/size(u[loaddofs],1)
time_storge[stp+1] = record[2]
if mod(stp,aa) == 0
n = Int(stp/aa)
numD[:,n].=d1[:]
numD2[:,:,n].=εᵖ_com
numD3[:,n].=u
numD4[:,:,n].=sigma
numD5[:,:,n].=εᵖ
numD6[:,:,n].= σc
numD7[:,:,n] .= Hn1
numD8[:,n] .= err_storage[:,n]
end
end
# STOnit[stp]=nit
# if nit==maxit && STOnit[stp-1] != maxit
# u_increment=0.25*u_increment
# elseif STOnit[stp]==maxit && STOnit[stp-1]==maxit && STOnit[stp-2]==maxit
# @info "Cannot converge"
# break
# end
end
return Fload1, Uload1, time_storge, numD, numD2, numD3, numD4, numD5, numD6, numD7, numD8
end
return monoli_initialStiff()
end