# Automatic difference of complex constitutive law

Hello, these days I am striving to realize a Jacobin of constitutive law in finite element simulation. I use the SparseDiffTools.jl, which accurately works at the linear elastic stage. However, when the displacements of nodes (finite element method) reaches the softening stage of the constitutive, I always encounter the error as below. I also marked the place where the error occurs in my function code. Any ideas are appreciated.

``````ERROR: UndefVarError: effectivestress not defined
Stacktrace:
[1] constitutive42(nodalforce1::Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(constitutive42), Float64}, Float64, 12}}, noddis::Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(constitutive42), Float64}, Float64, 12}})
@ Main c:\Users\DOJ14\OneDrive - University of Pittsburgh\PITT\Code\Julia\test_Mechanical\55
nonlinear LDPM_explicit working now 2\2\2.27 nonlinear LDPM_explicit working copy 4.jl:561
[2] forwarddiff_color_jacobian!(J::SparseArrays.SparseMatrixCSC{Float64, Int64}, f::typeof(constitutive42), x::Vector{Float64}, jac_cache::ForwardColorJacCache{Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(constitutive42), Float64}, Float64, 12}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(constitutive42), Float64}, Float64, 12}}, Vector{Float64}, Vector{Vector{NTuple{12, Float64}}}, Vector{Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}})
[3] forwarddiff_color_jacobian!(J::SparseArrays.SparseMatrixCSC{Float64, Int64}, f::Function, x::Vector{Float64}; dx::Vector{Float64}, colorvec::Vector{Int64}, sparsity::SparseArrays.SparseMatrixCSC{Float64, Int64})
[4] top-level scope
@ c:\Users\DOJ14\OneDrive - University of Pittsburgh\PITT\Code\Julia\test_Mechanical\55 nonlinear LDPM_explicit working now 2\2\2.27 nonlinear LDPM_explicit working copy 4.jl:1025
``````

My constitutive function is

``````function constitutive42(nodalforce1, noddis)

#get volumetric strain for every element
V_epsi = zeros(length1)
division = zeros(length1)

newverticecoor = verticeoftetrahedron1 + [noddis[(jj*6-5):(jj*6-3)] for jj = 1:length_nod][tttt1]
for tnumber = 1:length2
epsiVforte = (PhysicalMeshes.volume(PhysicalMeshes.Tetrahedron(PVector(newverticecoor[tnumber, 1]), PVector(newverticecoor[tnumber, 2]), PVector(newverticecoor[tnumber, 3]), PVector(newverticecoor[tnumber, 4]))) -
PhysicalMeshes.volume(PhysicalMeshes.Tetrahedron(PVector(verticeoftetrahedron1[tnumber, 1]), PVector(verticeoftetrahedron1[tnumber, 2]), PVector(verticeoftetrahedron1[tnumber, 3]), PVector(verticeoftetrahedron1[tnumber, 4])))) / PhysicalMeshes.volume(PhysicalMeshes.Tetrahedron(PVector(verticeoftetrahedron1[tnumber, 1]), PVector(verticeoftetrahedron1[tnumber, 2]), PVector(verticeoftetrahedron1[tnumber, 3]), PVector(verticeoftetrahedron1[tnumber, 4])))

for numele = 1:length1
if elements41[numele][1] ∈ tetrahedron1[tnumber, :] && elements41[numele][2] ∈ tetrahedron1[tnumber, :]
V_epsi = V_epsi + vcat(vec(zeros(numele - 1, 1)), epsiVforte, vec(zeros(length1 - numele, 1)))
division = division + vcat(vec(zeros(numele - 1, 1)), 1, vec(zeros(length1 - numele, 1)))
end
end
end
V_epsi_1 = deepcopy(@. V_epsi / division)
#get other strains for every element
for jii = 1:length_nod*6
nodalforce1[jii] = 0.0
end
for kjh in element_num_inter

epsiNk = BNjk2[kjh] * noddis[(6*elements41[kjh][2]-5):(6*elements41[kjh][2])] - BNik2[kjh] * noddis[(6*elements41[kjh][1]-5):(6*elements41[kjh][1])]
epsiMk = BMjk2[kjh] * noddis[(6*elements41[kjh][2]-5):(6*elements41[kjh][2])] - BMik2[kjh] * noddis[(6*elements41[kjh][1]-5):(6*elements41[kjh][1])]
epsiLk = BLjk2[kjh] * noddis[(6*elements41[kjh][2]-5):(6*elements41[kjh][2])] - BLik2[kjh] * noddis[(6*elements41[kjh][1]-5):(6*elements41[kjh][1])]
epsiTk = sqrt(epsiMk^2 + epsiLk^2)
##epsiNk
if epsiNk >= 0
effectiveepsi = sqrt(epsiTk^2 * alpha + epsiNk^2)
# tensile-shear boundary
omig = atan(epsiNk / sqrt(alpha) / epsiTk)
s = sin(omig)
c = cos(omig)
if omig == atan(sqrt(alpha) / mu)
effecstress0b = 0.5 * (sigmt + 2 * sigma) * sigmt / ((sigmt + sigma) * sin(omig))
else
effecstress0b = (-(sigmt + sigma) * s + sqrt(((sigmt + sigma) * s)^2 + (alpha * (c / mu)^2 - s^2) * (sigmt + 2 * sigma) * sigmt)) / (alpha * (c / mu)^2 - s^2)
end
Hs = 2 * alpha * EN / ((2 * alpha * EN * Gs) / (sigms^2 * le2[kjh]) - 1)
Ht = 2 * EN / ((2 * EN * Gt) / (sigmt^2 * le2[kjh]) - 1)
H = Hs + (Ht - Hs) * (2 * omig / pi)^nt
if maximum(effectivestrain_his[:, kjh]) == effectiveepsi
effectstressb = effecstress0b * exp(-H / effecstress0b * maximum([effectiveepsi - effecstress0b / EN, 0]))
else
effectstressb = effecstress0b * exp(-H / effecstress0b * maximum([maximum(effectivestrain_his[:, kjh]) - effecstress0b / EN, 0]))
end

if mode[t_step::Int64, kjh] == "continuous load"
effectivestress = minimum([effectiveepsi * EN, effectstressb])
effectivestress = minimum([maximum([effectivestress_unlo[t_step::Int64, kjh] - (effectivestrain_unlo[t_step::Int64, kjh] - effectiveepsi) * EN, 0+3*effectiveepsi]), effectstressb])
effectivestress = minimum([maximum([(effectiveepsi - kt * (effectivestrain_unlo[t_step::Int64, kjh] - effectivestress_unlo[t_step::Int64, kjh] / EN)) * EN, 0+3*effectiveepsi]), effectstressb])
end

#***!!!!This is exactly where the error happens!***
sigmNk = effectivestress / effectiveepsi * epsiNk

sigmTk = effectivestress / effectiveepsi * epsiTk * alpha

sigmLk = sigmTk * epsiLk / epsiTk

sigmMk = sigmTk * epsiMk / epsiTk

else #epsiNk < 0
#compressive

# rdv = (epsiNk - V_epsi_1[kjh]) / V_epsi_1[kjh]
# H = Hc / (1 + kapha2 * maximum([rdv-kapha1,0]))
# sigmcb = sigmc0 * exp(-H / sigmc0 * maximum([epsic0 - V_epsi_1[kjh], 0]))

# if mode_compre[t_step::Int64, kjh] == "continuous load"
rdv = (epsiNk - V_epsi_1[kjh]) / V_epsi_1[kjh]
H = Hc / (1 + kapha2 * maximum([rdv-kapha1,0]))
sigmcb = sigmc0 * exp(-H / sigmc0 * maximum([epsic0 - V_epsi_1[kjh], 0]))

if mode_compre[t_step::Int64, kjh] == "continuous load"
sigmNk = maximum([epsiNk * EN, sigmcb])
else
sigmNk = maximum([minimum([comprestress_unlo[t_step::Int64, kjh] - (comprestrain_unlo[t_step::Int64, kjh] - epsiNk) * Ed, sigmc0 - (comprestrain_unlo[t_step::Int64, kjh] - (comprestress_unlo[t_step::Int64, kjh] - sigmc0) / Ed - epsiNk) * EN, 0+3*epsiNk]), sigmcb])
end
# else
#     sigmNk = maximum([minimum([comprestress_unlo[t_step::Int64, kjh] - (comprestrain_unlo[t_step::Int64, kjh] - epsiNk) * Ed, sigmc0 - (comprestrain_unlo[t_step::Int64, kjh] - (comprestress_unlo[t_step::Int64, kjh] - sigmc0) / Ed - epsiNk) * EN, 0]), sigmcb])
# end##epsiNk = -0.0019340612033711502

#shear stress under compression

if sigmNk <= sigmN0
if sigmNk> sigmc0
sigmTb = sqrt((1-(sigmNk-sigmN0)^2/(sigmN0-sigmc0)^2)*sigmT0^2)
else
sigmTb =0.0
end
else
sigmTb = mu * sqrt((sigmNk - sigmt - sigma)^2 - sigma^2)
end

if t_step::Int64 == 1
sigmTk_f = epsiTk * ET
elseif fricstrain_his[t_step::Int64-1, kjh] == 300.0
sigmTk_f = shearsigm_his[t_step::Int64-1, kjh] + (epsiTk - sqrt(shearMstrain_his[t_step::Int64-1, kjh]^2 + shearLstrain_his[t_step::Int64-1, kjh]^2)) * ET
else
sigmTk_f = fricsigm_his[t_step::Int64-1, kjh] + (epsiTk - fricstrain_his[t_step::Int64-1, kjh]) * ET
end
if sigmTk_f > 0
sigmTk = minimum([sigmTk_f, sigmTb+3*epsiTk])
else
sigmTk = maximum([sigmTk_f, -sigmTb+3*epsiTk])
sigmTk = -sigmTk
end

if fricsigm_his[t_step::Int64-1, kjh] * fricsigm_his[t_step::Int64, kjh] < 0

epsiMk_1 = epsiMk - fricMstrain_his[t_step::Int64-1, kjh]# - (0 - fricsigm_his[t_step::Int64-1, kjh])/ET / sqrt(fricMstrain_his[t_step::Int64-1, kjh]^2 + fricLstrain_his[t_step::Int64-1, kjh]^2) * fricMstrain_his[t_step::Int64-1, kjh]
epsiLk_1 = epsiLk - fricLstrain_his[t_step::Int64-1, kjh]# - (0 - fricsigm_his[t_step::Int64-1, kjh])/ET / sqrt(fricMstrain_his[t_step::Int64-1, kjh]^2 + fricLstrain_his[t_step::Int64-1, kjh]^2) * fricLstrain_his[t_step::Int64-1, kjh]
else
epsiMk_1 = fricMstrain_reset_his[t_step::Int64-1, kjh] + epsiMk - fricMstrain_his[t_step::Int64-1, kjh]
epsiLk_1 = fricLstrain_reset_his[t_step::Int64-1, kjh] + epsiLk - fricLstrain_his[t_step::Int64-1, kjh]
end

sigmLk = sigmTk * epsiLk_1 / sqrt(epsiMk_1^2 + epsiLk_1^2) ##the ratio between L and M, use epsiLk / epsiTk maintains the relation between lk and mk

sigmMk = sigmTk * epsiMk_1 / sqrt(epsiMk_1^2 + epsiLk_1^2)

end
#println([sigmNk, sigmMk, sigmLk])
# if isnan(sigmNk) || isnan(sigmMk) || isnan(sigmLk)
#     println(kjh)
# end
# nodalforce1 = nodalforce1.*0.0 + vcat(vec(zeros(6*elements41[kjh][1]-6,1)), (-le2[kjh] * Ak2[kjh] * (sigmNk * BNik2[kjh] + sigmMk * BMik2[kjh] + sigmLk * BLik2[kjh]))',vec(zeros(length_nod*6 - (6*elements41[kjh][1]),1)))
# nodalforce1 = nodalforce1.*0.0 +[vec(zeros(1,6*elements41[kjh][2]-6))' (le2[kjh] * Ak2[kjh] * (sigmNk * BNjk2[kjh] + sigmMk * BMjk2[kjh] + sigmLk * BLjk2[kjh])) vec(zeros(1,length_nod*6 - (6*elements41[kjh][2])))']'
nodalforce1[(elements41[kjh][1]*6-5):(elements41[kjh][1]*6)] += (-le2[kjh] * Ak2[kjh] * (sigmNk * BNik2[kjh] + sigmMk * BMik2[kjh] + sigmLk * BLik2[kjh]))'
nodalforce1[(elements41[kjh][2]*6-5):(elements41[kjh][2]*6)] += (le2[kjh] * Ak2[kjh] * (sigmNk * BNjk2[kjh] + sigmMk * BMjk2[kjh] + sigmLk * BLjk2[kjh]))'

end
nothing

end
``````

I suspect none of the conditions in the if-else block matched. To test this, you can rewrite this as

``````effectivestress = if mode[t_step::Int64, kjh] == "continuous load"
minimum([effectiveepsi * EN, effectstressb])
minimum([maximum([effectivestress_unlo[t_step::Int64, kjh] - (effectivestrain_unlo[t_step::Int64, kjh] - effectiveepsi) * EN, 0+3*effectiveepsi]), effectstressb])
Then, if none of the conditions match `effectivestress` will be `nothing` and you can see what unexpected value is present.