fixed elasticityproblem bug

src/api_functions.jl

@@ -5,7 +5,6 @@ function add_boundary_condition!(case::Simulation, bc::NeumannBC)
     push!(case.neumann_boundary_conditions, bc)
 end
 
-
 """
 Add node to model and renumber for output
 """

src/api_types.jl

@@ -1,6 +1,7 @@
 # This file is a part of JuliaFEM.
 # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
 
+
 type NeumannBC
     set_name :: ASCIIString
     value :: Any

src/elasticity.jl

@@ -4,6 +4,7 @@
 # Elasticity problems
 
 abstract ElasticityProblem <: AbstractProblem
+abstract ElasticPlasticProblem <: AbstractProblem
 
 function get_unknown_field_name{P<:ElasticityProblem}(::Type{P})
     return "displacement"
@@ -13,8 +14,20 @@ function get_unknown_field_type{P<:ElasticityProblem}(::Type{P})
     return Vector{Float64}
 end
 
+function get_unknown_field_name{P<:ElasticPlasticProblem}(::Type{P})
+    return "displacement"
+end
+
+function get_unknown_field_type{P<:ElasticPlasticProblem}(::Type{P})
+    return Vector{Float64}
+end
+
 function ElasticityProblem(dim::Int=3, elements=[])
-    return Problem{PlaneStressElasticityProblem}(dim, elements)
+    return Problem{ElasticityProblem}(dim, elements)
+end
+
+function ElasticPlasticProblem(dim::Int=3, elements=[])
+    return Problem{ElasticPlasticProblem}(dim, elements)
 end
 
 abstract PlaneStressElasticityProblem <: ElasticityProblem
@@ -104,3 +117,60 @@ function get_residual_vector{P<:ElasticityProblem}(problem::Problem{P}, element:
     return vec(r)
 end
 
+"""
+
+"""
+function get_residual_vector{P<:ElasticPlasticProblem}(problem::Problem{P}, element::Element, ip::IntegrationPoint, time::Number; variation=nothing)
+
+
+    # u = element("displacement", ip, time, variation)
+
+    r = zeros(Float64, problem.dim, length(element))
+
+    # internal forces
+    if haskey(element, "youngs modulus") && haskey(element, "poissons ratio")
+        u = element("displacement", time, variation)
+        grad = element(ip, time, Val{:grad})
+#        gradu = element("displacement", ip, time, Val{:grad}, variation)
+        gradu = grad*u
+
+        F = I + gradu # deformation gradient
+
+        young = element("youngs modulus", ip, time)
+        poisson = element("poissons ratio", ip, time)
+
+        mu = young/(2*(1+poisson))
+        lambda = young*poisson/((1+poisson)*(1-2*poisson))
+        if P == PlaneStressElasticityProblem
+            lambda = 2*lambda*mu/(lambda + 2*mu)  # <- correction for 2d problems
+        end
+      
+
+        E = 1/2*(F'*F - I)       # large strain
+        #E = 1/2*(gradu + gradu')  # finite strain
+        S = lambda*trace(E)*I + 2*mu*E
+
+        #J = det(element, ip, time)
+        #T = J^-1*F*S*F'
+        #ip["cauchy stress"] = T
+        #ip["gl strain"] = E
+
+        r += F*S*grad
+    end
+
+    # external forces - volume load
+    if haskey(element, "displacement load")
+        basis = element(ip, time)
+        b = element("displacement load", ip, time)
+        r -= b*basis
+    end
+
+    # external forces - surface traction force
+    if haskey(element, "displacement traction force")
+        basis = element(ip, time)
+        T = element("displacement traction force", ip, time)
+        r -= T*basis
+    end
+
+    return vec(r)
+end

src/vonmises.jl

@@ -1,6 +1,7 @@
 using ForwardDiff
 using NLsolve
 
+
 function outer_prod(a, b)
     out = zeros(3,3,3,3)
     for i=1:3

test/test_von_mises_material.jl

@@ -13,20 +13,20 @@ function test_von_mises_basic()
     ν = 0.3
     C = stiffnessTensor(E, ν)
 
-    ϵ_tot = zeros(Float64, (steps, 6))
-    ϵ_tot2 = zeros(Float64, (steps, 6))
-    ϵ_tot3 = zeros(Float64, (steps, 6))
+    strain_tot = zeros(Float64, (steps, 6))
+    strain_tot2 = zeros(Float64, (steps, 6))
+    strain_tot3 = zeros(Float64, (steps, 6))
 
     # Adding only strain in x-axis and counting for the poisson effect
-    ϵ_tot[:, 1] = strain_max * sin(2 * pi * linspace(0, num_cycles, steps))
-    ϵ_tot[:, 2] = strain_max * sin(2 * pi * linspace(0, num_cycles, steps)).*-ν
-    ϵ_tot[:, 3] = strain_max * sin(2 * pi * linspace(0, num_cycles, steps)).*-ν
-    ϵ_tot[:, 4] = strain_max / 10 * sin(2 * pi * linspace(0, num_cycles, steps))
+    strain_tot[:, 1] = strain_max * sin(2 * pi * linspace(0, num_cycles, steps))
+    strain_tot[:, 2] = strain_max * sin(2 * pi * linspace(0, num_cycles, steps)).*-ν
+    strain_tot[:, 3] = strain_max * sin(2 * pi * linspace(0, num_cycles, steps)).*-ν
+    strain_tot[:, 4] = strain_max / 10 * sin(2 * pi * linspace(0, num_cycles, steps))
 
-    ϵ_last = zeros(Float64, (6))
-    ϵᵖ = zeros(Float64, (6))
-    σ = zeros(Float64, (6, 1))
-    σy =  200.0
+    strain_last = zeros(Float64, (6))
+    strain_p = zeros(Float64, (6))
+    stress = zeros(Float64, (6, 1))
+    stress_y =  200.0
     ss = Float64[]
     ee = Float64[]
 
@@ -47,19 +47,19 @@ function test_von_mises_basic()
         a[3, 2] = b[4]
     end
 
-    mat = State(C, σy, zeros(Float64, 6), zeros(Float64, 6))
+    mat = State(C, stress_y, zeros(Float64, 6), zeros(Float64, 6))
 
     info("Starting calculation")
     tic()
     for i=1:steps
-        ϵ_new = reshape(ϵ_tot[i, :, :], (6, 1))
-        dϵ = ϵ_new - mat.ϵ
-        calculate_stress!(dϵ, mat, Val{:vonMises})
-        mat.ϵ += vec(dϵ)
-        push!(ss, mat.σ[1])
-        push!(ee, mat.ϵ[1])
+        strain_new = reshape(strain_tot[i, :, :], (6, 1))
+        dstrain = strain_new - mat.strain
+        calculate_stress!(dstrain, mat, Val{:vonMises})
+        mat.strain += vec(dstrain)
+        push!(ss, mat.stress[1])
+        push!(ee, mat.strain[1])
 
-        fill_tensor(eig_stress, mat.σ)
+        fill_tensor(eig_stress, mat.stress)
         eig_vals[i, :] = sort(eigvals(eig_stress))
     end
     toc()