diff --git a/src/vonMises.jl b/src/vonMises.jl index 405b8ea..f72def5 100644 --- a/src/vonMises.jl +++ b/src/vonMises.jl @@ -24,22 +24,56 @@ Returns Array{Float64, (6,6)} """ #= -my_kron(i,j) = i == j ? 1 : 0 -II = zeros(Float64, (3, 3, 3, 3)) -for i=1:3 - for j=1:3 - for k=1:3 - for l=1:3 - A[i,j,k,l] = my_kron(i, j) * my_kron(k, l) +function outer_prod(a, b) + out = zeros(3,3,3,3) + for i=1:3 + for j=1:3 + for k=1:3 + for l=1:3 + out[i, j, k, l] = a[i, j] * b[k, l] + end end end end + out +end + +function double_contract(a, b) + out = zeros(3, 3) + for i=1:3 + for j=1:3 + for k=1:3 + for l=1:3 + out[i, j] = a[i,j,k,l] * b[k, l] + end + end + end + end + out +end + +function stiffnessTensor(youngs_modulus, poissons_ratio, ::Type{Val{:isotropic}}) + E = youngs_modulus + v = poissons_ratio + my_kron(i,j) = i == j ? 1 : 0 + II = zeros(Float64, (3, 3, 3, 3)) + for i=1:3 + for j=1:3 + for k=1:3 + for l=1:3 + II[i,j,k,l] = my_kron(i, k) * my_kron(j, l) + end + end + end + end + mu = E/(2*(1+v)) + lambda = E*v/((1+v)*(1-2*v)) + I = eye(3) + return lambda * outer_prod(I, I) + 2 * mu * II + #C = λ * I ⊗ I + 2 * μ * II end =# -#function outer_product(a, b) -# -#end -function hookeStiffnessTensor(E, ν) +function stiffnessTensor(E, ν) a = 1 - ν b = 1 - 2*ν c = 1 + ν @@ -52,14 +86,6 @@ function hookeStiffnessTensor(E, ν) 0 0 0 0 0 b].*multiplier end -# Pick material values -#E = 200.0e3 -#mu = 0.3 -#C = hookeStiffnessTensor(E, ν) -#lambda = -#nu = -#C = λ * I ⊗ I + 2 * μ * II - type State C :: Array{Float64, 2} σ_y :: Float64 diff --git a/test/test_von_mises_material.jl b/test/test_von_mises_material.jl index 5829a81..45de507 100644 --- a/test/test_von_mises_material.jl +++ b/test/test_von_mises_material.jl @@ -1,13 +1,7 @@ module VonMisesTests using PyPlot - -macro R_str(s) - s -end - - -using JuliaFEM.MaterialModels: hookeStiffnessTensor, calculate_stress!, State +using JuliaFEM.MaterialModels: stiffnessTensor, calculate_stress!, State function test_von_mises_basic() @@ -17,7 +11,7 @@ function test_von_mises_basic() E = 200.0e3 nu = 0.3 ν = 0.3 - C = hookeStiffnessTensor(E, ν) + C = stiffnessTensor(E, ν) ϵ_tot = zeros(Float64, (steps, 6)) ϵ_tot2 = zeros(Float64, (steps, 6))