Closed
Description
The following code gives different results for the LinearAlgebra det and a hand-written determinant function:
MWE
using LinearAlgebra
using ForwardDiff: derivative, gradient, jacobian, hessian
E, nu = (7.0e6, 0.3)
lambda = E * nu / (1 + nu) / (1 - 2*(nu));
mu = E / 2. / (1+nu);
function strain6vto3x3t!(t::Array{T,2}, v::Array{T,1}) where {T}
t[1,1] = v[1];
t[2,2] = v[2];
t[3,3] = v[3];
t[1,2] = v[4]/2.;
t[2,1] = v[4]/2.;
t[1,3] = v[5]/2.;
t[3,1] = v[5]/2.;
t[3,2] = v[6]/2.;
t[2,3] = v[6]/2.;
return t
end
function strain3x3tto6v!(v::Array{T,1}, t::Array{T,2}) where {T}
v[1] = t[1,1];
v[2] = t[2,2];
v[3] = t[3,3];
v[4] = t[1,2] + t[2,1];
v[5] = t[1,3] + t[3,1];
v[6] = t[3,2] + t[2,3];
return v
end
mydet(C) = begin
C[1,1] * C[2,2] * C[3,3] +
C[1,2] * C[2,3] * C[3,1] +
C[1,3] * C[2,1] * C[3,2] -
C[1,3] * C[2,2] * C[3,1] -
C[1,2] * C[2,1] * C[3,3] -
C[1,1] * C[2,3] * C[3,2]
end
for i in 1:100
a = rand(3, 3)
@assert (det(a) - mydet(a)) / det(a) <= 1.0e-6
end
function strainenergy(Cv, lambda, mu)
C = fill(zero(eltype(Cv)), 3, 3)
strain6vto3x3t!(C, Cv)
J = sqrt(mydet(C))
lJ = log(J)
return mu/2*(tr(C)-3) - mu*lJ + (lambda/2)*lJ^2;
end
function strainenergy1(Cv, lambda, mu)
C = fill(zero(eltype(Cv)), 3, 3)
strain6vto3x3t!(C, Cv)
J = sqrt(det(C))
lJ = log(J)
return mu/2*(tr(C)-3) - mu*lJ + (lambda/2)*lJ^2;
end
Fn1 = [1.1 0 0; 0 1.0 0; 0 0 1.0]
Fn = [1.0 0 0; 0 1.0 0; 0 0 1.0]
C = Fn1'*Fn1;
Cv = fill(0.0, 6)
strain3x3tto6v!(Cv, C)
@show Da = hessian(Cv -> strainenergy(Cv, lambda, mu), Cv);
@show Da = hessian(Cv -> strainenergy1(Cv, lambda, mu), Cv);
Julia Version 1.3.0-alpha.0
Commit 6c11e7c2c4 (2019-07-23 01:46 UTC)
Platform Info:
OS: Windows (x86_64-w64-mingw32)
CPU: Intel(R) Core(TM) i7-8705G CPU @ 3.10GHz
WORD_SIZE: 64
LIBM: libopenlibm
LLVM: libLLVM-6.0.1 (ORCJIT, skylake)
Status `C:\Users\PetrKrysl\.julia\environments\v1.3\Project.toml`
[f6369f11] ForwardDiff v0.10.3
Metadata
Metadata
Assignees
Labels
No labels