I am working on a modified finite element simulation. I have a complicated function to project a 648 units long vector to a 648 units long vector (i.e., from nodal displacements to nodal forces, the magnitude will be maginified to 10000-to-10000). The function (project) is kind of complicated, I put it as follows. My goal is to get the Jacobian matirix of the Vector-to-Vector function at every Vector input. I have tried some automatic differentiation packages, as well as Zygote, but failed always. Every effective recommendation or clue is appreciated.
I also upload the data (the constants in the fucntion) needed as a jld2 in my Github GitHub - DonggeJia/the-function-needs-Jacobian: A project from nodal displacement to nodal force of LDPM
using PhysicalMeshes
function constitutive(noddis)
#get volumetric strain for every element
V_epsi = vec(zeros(length1, 1))
division = vec(zeros(length1, 1))
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
nodalforce1 = vec(zeros(length(pointsfinal)*6,1))
for kjh = 1:length1
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)
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) / (sigmt^2 * le2[kjh]) - 1)
Ht = 2 * EN / ((2 * EN * Gt) / (sigmt^2 * le2[kjh]) - 1)
H = Hs + (Ht - Hs) * (2 * omig / pi)^nt
effectstressb = effecstress0b * exp(-H / effecstress0b * maximum([maxeffecsigm[kjh] - effecstress0b / EN, 0]))
if effectiveepsi * EN <= effectstressb
effectivestress = effectiveepsi * EN
else
effectivestress = effectstressb
end
if epsiNk == 0.0
sigmNk = 0.0
else
sigmNk = effectivestress / effectiveepsi * epsiNk
end
sigmTk = effectivestress / effectiveepsi * epsiTk * alpha
if epsiLk == 0.0
sigmLk = 0.0
else
sigmLk = sigmTk * epsiLk / epsiTk ##the ratio between L and M, use epsiLk / epsiTk maintains the relation between lk and mk
end
if epsiMk == 0.0
sigmMk = 0.0
else
sigmMk = sigmTk * epsiMk / epsiTk
end
else #epsiNk < 0
#compressive
rdv = (epsiNk - V_epsi_1[kjh]) / V_epsi_1[kjh]
if V_epsi_1[kjh] >= 0
if epsiNk >= sigmc0 / EN
sigmNk = epsiNk * EN
else
sigmNk = sigmc0
end
else #V_epsi_1[kjh] < 0
H = Hc / (1 + kapha * rdv)
sigmcb = sigmc0 * exp(H / sigmc0 * maximum([epsic0 - V_epsi_1[kjh], 0]))
if epsiNk >= sigmcb / EN
sigmNk = epsiNk * EN
else
sigmNk = sigmcb
end
end
#shear stress under compression
if sigmNk <= sigmN0
omig = atan(epsiNk / sqrt(alpha) / epsiTk)
s = sin(omig)
c = cos(omig)
sigmTb0 = -sigmc0 / sqrt(s^2 + alpha * c^2 / betta)
else
sigmTb0 = mu * sqrt((epsiNk - sigmt - sigma)^2 - sigma^2)
end
epsiT0 = sigmTb0 / ET
omig = atan(epsiNk / sqrt(alpha) / epsiTk)
Hs = 2 * alpha * EN / ((2 * alpha * EN * Gs) / (sigmt^2 * le2[kjh]) - 1)
H = Hs - Hs * (-2 * omig / pi)^nc
sigmTb = sigmTb0 * exp(-H / sigmTb0 * maximum([epsiTk - epsiT0, 0]))
if epsiTk < sigmTb / ET
sigmTk = epsiTk * ET
else
sigmTk = sigmTb
end
if epsiLk == 0.0
sigmLk = 0.0
else
sigmLk = sigmTk * epsiLk / epsiTk ##the ratio between L and M, use epsiLk / epsiTk maintains the relation between lk and mk
end
if epsiMk == 0.0
sigmMk = 0.0
else
sigmMk = sigmTk * epsiMk / epsiTk
end
end
nodalforce1 = nodalforce1 + 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 +[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])))']'
end
return nodalforce1
end