diff --git a/src/elements.jl b/src/elements.jl index 4183c81..5925559 100644 --- a/src/elements.jl +++ b/src/elements.jl @@ -9,15 +9,14 @@ type Element{E<:AbstractElement} integration_points :: Vector{IP} fields :: Dict{AbstractString, Field} properties :: E - dev :: Dict{Any, Any} end function Element{E<:AbstractElement}(::Type{E}, id::Int64, connectivity=[]) - return Element{E}(id, connectivity, [], Dict(), E(), Dict{Any, Any}()) + return Element{E}(id, connectivity, [], Dict(), E()) end function Element{E<:AbstractElement}(::Type{E}, connectivity=[]) - return Element{E}(-1, connectivity, [], Dict(), E(), Dict{Any, Any}()) + return Element{E}(-1, connectivity, [], Dict(), E()) end function getindex(element::Element, field_name::AbstractString) diff --git a/src/materials_plasticity.jl b/src/materials_plasticity.jl index 128165c..d1c000a 100644 --- a/src/materials_plasticity.jl +++ b/src/materials_plasticity.jl @@ -72,7 +72,7 @@ function radial_return(params, dstrain, D, stress_y, stress_base, yield_surface_ [vec(function_1); function_2] end -function ideal_plasticity!(stress_new, stress_last, dstrain_vec, D, params, Dtan, yield_surface_, type_) +function ideal_plasticity!(stress_new, stress_last, dstrain_vec, D, params, Dtan, yield_surface_, time, dt, type_) # Test stress dstress = vec(D * dstrain_vec) stress_trial = stress_last + dstress @@ -115,7 +115,5 @@ function ideal_plasticity!(stress_new, stress_last, dstrain_vec, D, params, Dtan dfds = dfds_(stress_new) Dtan[:,:] = Dc - (Dc * dfds * dfds' * Dc) / (dfds' * Dc * dfds)[1] - println("equivalent stress") - println(equivalent_stress(stress_new, type_)) end end diff --git a/src/problems_elasticity.jl b/src/problems_elasticity.jl index c4aa901..764125a 100644 --- a/src/problems_elasticity.jl +++ b/src/problems_elasticity.jl @@ -71,22 +71,34 @@ typealias Elasticity2DVolumeElements Union{Tri3, Tri6, Quad4, Quad8, Quad9} typealias Elasticity3DSurfaceElements Union{Poi1, Tri3, Tri6, Quad4, Quad8, Quad9} typealias Elasticity3DVolumeElements Union{Tet4, Wedge6, Hex8, Tet10, Hex20, Hex27} -function get_internal_params(params, ip_id, ::Type{Val{:type_2d}}) - if !(ip_id in keys(params)) - params[ip_id] = Dict{Any, Any}() - params[ip_id]["last_stress"] = [0.0,0.0,0.0] - params[ip_id]["last_strain"] = [0.0,0.0,0.0] +function initialize_internal_params!(params, ip, ::Type{Val{:type_2d}}) + param_keys = keys(params) + all_keys = ip.fields.keys + ip_fields = filter(x->isdefined(all_keys, x), collect(1:length(all_keys))) + + if !("params_initialized" in ip_fields) + for key in param_keys + update!(ip, key, 0.0 => params[key]) + end + update!(ip, "stress", 0.0 => [0.0,0.0,0.0]) + update!(ip, "strain", 0.0 => [0.0,0.0,0.0]) + update!(ip, "prev_time", 0.0 => 0.0) + update!(ip, "params_initialized", 0.0 => true) end - return (params[ip_id]["last_stress"], params[ip_id]["last_strain"]) end -function get_internal_params(params, ip_id, ::Type{Val{:type_3d}}) +function get_keys(element) + all_keys = element.fields.keys + idx = filter(x->isdefined(all_keys, x), collect(1:length(all_keys))) + map(x -> all_keys[x], idx) +end + +function initialize_internal_params!(params, ip_id, ::Type{Val{:type_3d}}) if !(ip_id in keys(params)) params[ip_id] = Dict{Any, Any}() params[ip_id]["last_stress"] = [0.0,0.0,0.0,0.0,0.0,0.0] params[ip_id]["last_strain"] = [0.0,0.0,0.0,0.0,0.0,0.0] end - return (params[ip_id]["last_stress"], params[ip_id]["last_strain"]) end """ Elasticity equations for 2d cases. """ @@ -153,16 +165,34 @@ function assemble{El<:Elasticity2DVolumeElements}(problem::Problem{Elasticity}, else error("unknown plane formulation: $(props.formulation)") end + # calculate stress - if "plasticity" in keys(element.dev) - plastic_def = element.dev["plasticity"] - calculate_stress! = plastic_def["stress"] + element_keys = get_keys(element) + + if "plasticity" in element_keys + plastic_def = element("plasticity")[ip.id] + + calculate_stress! = plastic_def["type"] yield_surface_ = plastic_def["yield_surface"] params = plastic_def["params"] - (stress_last, strain_last) = get_internal_params(element.dev, ip.id, Val{:type_2d}) + + initialize_internal_params!(params, ip, Val{:type_2d}) + + if time == 0.0 + error("Given step time = $(time). Please select time > 0.0") + end + + t_last = ip("prev_time", time) + update!(ip, "prev_time", time => t_last) + + dt = time - t_last + + stress_last = ip("stress", t_last) + strain_last = ip("strain", t_last) + dstrain_vec = strain_vec - strain_last stress_vec = [0.0, 0.0, 0.0] - calculate_stress!(stress_vec, stress_last, dstrain_vec, D, params, Dtan, yield_surface_, Val{:type_2d}) + calculate_stress!(stress_vec, stress_last, dstrain_vec, D, params, Dtan, yield_surface_, time, dt, Val{:type_2d}) else stress_vec = D * ([1.0, 1.0, 2.0] .* strain_vec) Dtan[:,:] = D[:,:] @@ -176,6 +206,7 @@ function assemble{El<:Elasticity2DVolumeElements}(problem::Problem{Elasticity}, Km += w*BL'*Dtan*BL + # stress = [stress_vec[1] stress_vec[3]; stress_vec[3] stress_vec[2]] # cauchy_stress = F'*stress*F/det(F) # cauchy_stress = [cauchy_stress[1,1]; cauchy_stress[2,2]; cauchy_stress[1,2]] diff --git a/test/test_elasticplastic_2d_nonhomogenious_boundary_conditions.jl b/test/test_elasticplastic_2d_nonhomogenious_boundary_conditions.jl index 60c6ba5..5eee961 100644 --- a/test/test_elasticplastic_2d_nonhomogenious_boundary_conditions.jl +++ b/test/test_elasticplastic_2d_nonhomogenious_boundary_conditions.jl @@ -23,9 +23,13 @@ using JuliaFEM.Testing update!(element, "geometry", nodes) update!(element, "youngs modulus", 288.0) update!(element, "poissons ratio", 1/3) - element.dev["plasticity"] = Dict{Any, Any}("stress" => JuliaFEM.ideal_plasticity!, - "yield_surface" => Val{:von_mises}, - "params" => Dict("yield_stress" => 175.0)) + + plastic_parameters = Dict{Any, Any}("type" => JuliaFEM.ideal_plasticity!, + "yield_surface" => Val{:von_mises}, + "params" => Dict("yield_stress" => 175.0)) + to_integ_points = Dict() + map(x-> to_integ_points[x] = plastic_parameters, get_connectivity(element)) + update!(element, "plasticity", to_integ_points) push!(block, element) # boundary conditions @@ -40,6 +44,7 @@ using JuliaFEM.Testing push!(bc, bel1, bel2, bel3) solver = NonlinearSolver("solve block problem") + solver.time = 1.0 push!(solver, block, bc) solver()