diff --git a/src/JuliaFEM.jl b/src/JuliaFEM.jl
index 3240258..f5a25f0 100644
--- a/src/JuliaFEM.jl
+++ b/src/JuliaFEM.jl
@@ -67,6 +67,9 @@ export Problem, AbstractProblem, FieldProblem, BoundaryProblem,
 include("problems_elasticity.jl")
 export Elasticity
 
+include("materials_plasticity.jl")
+export plastic_von_mises
+
 include("problems_dirichlet.jl")
 export Dirichlet
 
diff --git a/src/elements.jl b/src/elements.jl
index 59b9c0b..4183c81 100644
--- a/src/elements.jl
+++ b/src/elements.jl
@@ -9,14 +9,15 @@ 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())
+    return Element{E}(id, connectivity, [], Dict(), E(), Dict{Any, Any}())
 end
 
 function Element{E<:AbstractElement}(::Type{E}, connectivity=[])
-    return Element{E}(-1, connectivity, [], Dict(), E())
+    return Element{E}(-1, connectivity, [], Dict(), E(), Dict{Any, Any}())
 end
 
 function getindex(element::Element, field_name::AbstractString)
@@ -71,7 +72,7 @@ julia> el([0.0, 0.0], 0.0, 1)
 
 julia> el([0.0, 0.0], 0.0, 2)
 2x8 Array{Float64,2}:
- 0.25  0.0   0.25  0.0   0.25  0.0   0.25  0.0 
+ 0.25  0.0   0.25  0.0   0.25  0.0   0.25  0.0
  0.0   0.25  0.0   0.25  0.0   0.25  0.0   0.25
 
 """
diff --git a/src/vonmises.jl b/src/materials_plasticity.jl
similarity index 76%
rename from src/vonmises.jl
rename to src/materials_plasticity.jl
index ecc477f..257737d 100644
--- a/src/vonmises.jl
+++ b/src/materials_plasticity.jl
@@ -1,45 +1,18 @@
 using ForwardDiff
-
-"""
-Create a isotropic Hooke material matrix C
-
-More information: http://www.efunda.com/formulae/solid_mechanics/mat_mechanics/hooke_isotropic.cfm
-                  https://en.wikipedia.org/wiki/Hooke's_law
-                  http://www.ce.berkeley.edu/~sanjay/ce231mse211/symidentity.pdf
-Parameters
-----------
-    E: Float
-        Elastic modulus
-    ν: Float
-        Poisson constant
-
-Returns
--------
-    Array{Float64, (6,6)}
-"""
-function stiffnessTensor(E, ν)
-    a = 1 - ν
-    b = 1 - 2*ν
-    c = 1 + ν
-    multiplier = E / (b * c)
-    return Float64[a ν ν 0 0 0;
-                   ν a ν 0 0 0;
-                   ν ν a 0 0 0;
-                   0 0 0 b 0 0;
-                   0 0 0 0 b 0;
-                   0 0 0 0 0 b].*multiplier
-end
+using NLsolve
 
 # Creating functions for newton: xₙ₊₁ = xₙ - df⁻¹ * f and initial values
-function find_root!(f, df, x; max_iter=50, norm_acc=1e-10)
+function find_root!(f, df, x; max_iter=50, norm_acc=1e-9)
     converged = false
     for i=1:max_iter
-        dx = df(x) \ -f(x)
-        x += dx
+        dx = -df(x) \ f(x)
         norm(dx) < norm_acc && (converged = true; break)
+
+        x += dx
+
     end
     converged || error("no convergence!")
-    x
+    return x
 end
 
 type State
@@ -227,15 +200,6 @@ end
 """
 http://www.efunda.com/formulae/solid_mechanics/mat_mechanics/hooke_plane_stress.cfm
 """
-function stiffnessTensorPlaneStress(E, ν)
-    a = 1 - ν^2
-    b = 1 - ν
-    multiplier = E / a
-    return Float64[1 ν 0;
-                   ν 1 0;
-                   0 0 b].*multiplier
-end
-
 # von mises: plane stress
 # https://andriandriyana.files.wordpress.com/2008/03/yield_criteria.pdf
 function stress_eq_plane_stress(stress)
@@ -244,6 +208,7 @@ function stress_eq_plane_stress(stress)
     # http://www.engineersedge.com/material_science/principal_vonmises_stress__13418.htm
     se1 = (s1 + s2)/2 + sqrt(((s1 - s2)/2)^2 + t12^2)
     se2 = (s1 + s2)/2 - sqrt(((s1 - s2)/2)^2 + t12^2)
+
     return sqrt(se1^2 -se1*se2 + se2^2)
 end
 
@@ -252,53 +217,65 @@ function vonMisesYieldPlaneStress(stress, stress_y)
     stress_eq_plane_stress(stress) - stress_y
 end
 
-function vonMisesRootPlaneStress(params, dstrain, C, stress_y, stress_base)
+function vonMisesRootPlaneStress(params, dstrain, D, stress_y, stress_base)
 
     # Creating wrapper for gradient
     vm_wrap(stress_) = vonMisesYieldPlaneStress(stress_, stress_y)
-    dfds = ForwardDiff.gradient(vm_wrap)
+    dfds = x -> ForwardDiff.gradient(vm_wrap, x)
 
     # Stress rate and total strain
     dstress = params[1:3]
-    stress_tot = vec(stress_base) + params[1:3]
+    stress_tot = stress_base + params[1:3]
 
     # Calculating plastic strain rate
     dstrain_p = params[end] * dfds(stress_tot)
 
     # Calculating equations
-    function_1 = dstress - C * (dstrain - dstrain_p)
+    function_1 = dstress - D * (dstrain - dstrain_p)
     function_2 = vm_wrap(stress_tot)
     [vec(function_1); function_2]
 end
 
-function calculate_stress(dstrain, stress, C, stress_y,
-                          ::Type{Val{:vonMises}},
-                          ::Type{Val{:PlaneStressElasticPlasticProblem}})
+# http://homes.civil.aau.dk/lda/continuum/plast.pdf
+function plastic_von_mises!(stress, dstrain_vec, D, params, Dtan)
     # Test stress
-    dstress = C * dstrain
+    dstress = vec(D * dstrain_vec)
     stress_tria = stress + dstress
+    stress_y = params["yield_stress"]
 
     # Calculating and checking for yield
     yield = vonMisesYieldPlaneStress(stress_tria, stress_y)
     if isless(yield, 0.0)
-        return dstress, zeros(3)
+        stress[:] = stress_tria[:]
+        Dtan[:,:] = D[:,:]
     else
-        info("yielded")
-        # Yielding happened
-        # Creating functions for newton: xₙ₊₁ = xₙ - df⁻¹ * f and initial values
+        # Creating functions for newton: xₙ₊₁ = xₙ - df⁻¹ \ f and initial values
         x = [vec(stress_tria - stress); 0.0]
-        f(stress_)  = vonMisesRootPlaneStress(stress_, dstrain, C, stress_y, stress)
-        df = ForwardDiff.jacobian(f)
+        f = stress_ -> vonMisesRootPlaneStress(stress_, dstrain_vec, D, stress_y, stress)
+        df = x -> ForwardDiff.jacobian(f, x)
+
         # Calculating root
-        results = find_root!(f, df, x)
+        results = nlsolve(not_in_place(f), x).zero
         dstress = results[1:3]
-        stress_tot = stress + dstress
+
+        stress_new = stress + dstress
         plastic_multiplier = results[end]
-        vm_wrap(stress_) = vonMisesYieldPlaneStress(stress_, stress_y)
-        dfds = ForwardDiff.gradient(vm_wrap)
-        dep = plastic_multiplier * dfds(vec(stress_tot))
-        info("II ", stress_tot)
-        info(vm_wrap(stress_tot))
-        return dstress, dep
+        f_ = stress_ -> vonMisesYieldPlaneStress(stress_, stress_y)
+        dfds_ = x -> ForwardDiff.gradient(f_, x)
+        dep = plastic_multiplier * dfds_(vec(stress_new))
+
+        D2g = x -> ForwardDiff.hessian(f_, x)
+        Dc = (D^-1 + plastic_multiplier * D2g(stress_new))^-1
+        dfds = dfds_(stress_new)
+        Dtan = Dc - (Dc * dfds * dfds' * Dc) / (dfds' * Dc * dfds)[1]
+        println("plastic stress")
+        println(stress_new)
+        println(Dtan)
+        stress[:] = stress_new[:]
+        # stress[:] = D * dstrain_vec
+        println("elastic stress")
+        println(stress)
+        println(D)
+        Dtan[:,:] = D[:,:]
     end
 end
diff --git a/src/problems_elasticity.jl b/src/problems_elasticity.jl
index b7dc871..2aad6ab 100644
--- a/src/problems_elasticity.jl
+++ b/src/problems_elasticity.jl
@@ -71,6 +71,14 @@ 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{:planestress}})
+    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]
+    end
+    return (params[ip_id]["last_stress"], params[ip_id]["last_strain"])
+end
 
 """ Elasticity equations for 2d cases. """
 function assemble{El<:Elasticity2DVolumeElements}(problem::Problem{Elasticity}, element::Element{El}, time, ::Type{Val{:plane}})
@@ -137,7 +145,22 @@ function assemble{El<:Elasticity2DVolumeElements}(problem::Problem{Elasticity},
             error("unknown plane formulation: $(props.formulation)")
         end
         # calculate stress
-        stress_vec = D * ([1.0, 1.0, 2.0] .* strain_vec)
+        if "plasticity" in keys(element.dev)
+            plastic_def = element.dev["plasticity"]
+            calculate_stress! = plastic_def["stress"]
+            params = plastic_def["params"]
+            (stress_last, strain_last) = get_internal_params(element.dev, ip.id, Val{:planestress})
+            dstrain_vec = strain_vec - strain_last
+            Dtan = [0.0 0.0 0.0;
+                    0.0 0.0 0.0;
+                    0.0 0.0 0.0]
+
+            calculate_stress!(stress_last, dstrain_vec, D, params, Dtan)
+            stress_vec = stress_last
+        else
+            stress_vec = D * ([1.0, 1.0, 2.0] .* strain_vec)
+            Dtan = D
+        end
 
         :strain in props.store_fields && update!(ip, "strain", time => strain_vec)
         :stress in props.store_fields && update!(ip, "stress", time => stress_vec)
@@ -145,7 +168,7 @@ function assemble{El<:Elasticity2DVolumeElements}(problem::Problem{Elasticity},
         :stress22 in props.store_fields && update!(ip, "stress22", time => stress_vec[2])
         :stress12 in props.store_fields && update!(ip, "stress12", time => stress_vec[3])
 
-        Km += w*BL'*D*BL
+        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)
diff --git a/test/test_elasticity_2d_nonhomogeneous_boundary_conditions.jl b/test/test_elasticity_2d_nonhomogeneous_boundary_conditions.jl
index d4d3754..268fe87 100644
--- a/test/test_elasticity_2d_nonhomogeneous_boundary_conditions.jl
+++ b/test/test_elasticity_2d_nonhomogeneous_boundary_conditions.jl
@@ -70,4 +70,3 @@ using JuliaFEM.Testing
     end
 =#
 end
-
diff --git a/test/test_elasticplastic_2d_nonhomogenious_boundary_conditions.jl b/test/test_elasticplastic_2d_nonhomogenious_boundary_conditions.jl
new file mode 100644
index 0000000..f8a1260
--- /dev/null
+++ b/test/test_elasticplastic_2d_nonhomogenious_boundary_conditions.jl
@@ -0,0 +1,53 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
+
+using JuliaFEM
+using JuliaFEM.Preprocess
+using JuliaFEM.Testing
+
+# @testset "2d nonlinear elasticity: test nonhomogeneous boundary conditions and stress calculation" begin
+
+    # field problem
+    block = Problem(Elasticity, "BLOCK", 2)
+    block.properties.formulation = :plane_stress
+    block.properties.finite_strain = true
+    block.properties.geometric_stiffness = true
+
+    nodes = Dict{Int, Vector{Float64}}(
+        1 => [0.0, 0.0],
+        2 => [1.0, 0.0],
+        3 => [1.0, 1.0],
+        4 => [0.0, 1.0])
+
+    element = Element(Quad4, [1, 2, 3, 4])
+    update!(element, "geometry", nodes)
+    update!(element, "youngs modulus", 288.0)
+    update!(element, "poissons ratio", 1/3)
+    element.dev["plasticity"] = Dict{Any, Any}("stress" => JuliaFEM.plastic_von_mises!,
+                                               "params" => Dict("yield_stress" => 175.0))
+    push!(block, element)
+
+    # boundary conditions
+    bc = Problem(Dirichlet, "bc", 2, "displacement")
+    bel1 = Element(Seg2, [1, 2])
+    bel2 = Element(Seg2, [3, 4])
+    bel3 = Element(Seg2, [4, 1])
+    update!([bel1, bel2, bel3], "geometry", nodes)
+    update!(bel1, "displacement 2", 0.0)
+    update!(bel2, "displacement 2", 0.5)
+    update!(bel3, "displacement 1", 0.0)
+    push!(bc, bel1, bel2, bel3)
+
+    solver = NonlinearSolver("solve block problem")
+    push!(solver, block, bc)
+    solver()
+
+    # from code aster
+    eps_expected = [-2.08333312468287E-01, 6.25000000000000E-01, 0.0]
+    sig_expected = [ 4.50685020821470E-06, 4.62857140373777E+02, 0.0]
+    u3_expected =  [-2.36237356855269E-01, 5.00000000000000E-01]
+
+    u3 = reshape(block.assembly.u, 2, 4)[:, 3]
+    info("u3 = $u3")
+    #@test isapprox(u3, u3_expected, atol=1.0e-5)
+# end
diff --git a/test/test_von_mises_material.jl b/test/test_von_mises_material.jl
index 27012b9..601401b 100644
--- a/test/test_von_mises_material.jl
+++ b/test/test_von_mises_material.jl
@@ -1,145 +1,149 @@
 
-#using PyPlot
+using PyPlot
+using JuliaFEM
 using JuliaFEM.Testing
 #using JuliaFEM.MaterialModels: stiffnessTensor, calculate_stress, State
 #using JuliaFEM.MaterialModels: stiffnessTensorPlaneStress
 
 
-function test_von_mises_3D_basic()
+# function test_von_mises_3D_basic()
+#
+#     steps = 1000
+#     strain_max = 0.003
+#     num_cycles = 3
+#     E = 200.0e3
+#     nu =  0.3
+#     ν = 0.3
+#     C = stiffnessTensor(E, ν)
+#
+#     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
+#     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))
+#
+#     strain_last = zeros(Float64, (6))
+#     strain_p = zeros(Float64, (6))
+#     stress = zeros(Float64, (6, 1))
+#     stress_y =  200.0
+#     ss = Float64[]
+#     ee = Float64[]
+#
+#     eig_stress = zeros(Float64, (3, 3))
+#     eig_vals = zeros(Float64, (steps, 3))
+#
+#     function fill_tensor(a, b)
+#         a[1, 1] = b[1]
+#         a[2, 2] = b[2]
+#         a[3, 3] = b[3]
+#
+#         a[1, 2] = b[6]
+#         a[1, 3] = b[5]
+#         a[2, 3] = b[4]
+#
+#         a[2, 1] = b[6]
+#         a[3, 1] = b[5]
+#         a[3, 2] = b[4]
+#     end
+#
+#     mat = State(C, stress_y, zeros(Float64, 6), zeros(Float64, 6))
+#
+#     info("Starting calculation")
+#     tic()
+#     #=
+#     for i=1:steps
+#         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.stress)
+#         eig_vals[i, :] = sort(eigvals(eig_stress))
+#     end
+#     =#
+#     stress = zeros(Float64, 6)
+#     strain = zeros(Float64, 6)
+#     for i=1:steps
+#         strain_new = reshape(strain_tot[i, :, :], (6, 1))
+#         dstrain = strain_new - strain
+#         calculate_stress!(dstrain, stress, C, stress_y, Val{:vonMises})
+#         strain = vec(strain_new)
+#         push!(ss, stress[1])
+#         push!(ee, strain[1])
+#         fill_tensor(eig_stress, stress)
+#         eig_vals[i, :] = sort(eigvals(eig_stress))
+#     end
+#
+#     toc()
+#     # ================ Plotting =================== #
+#     n(θ, ϕ) = [sin(θ)*cos(ϕ)
+#                sin(θ)*sin(ϕ)
+#                cos(θ)]
+#     m(θ, ϕ, χ) = [-sin(ϕ)*cos(χ)-cos(θ)*cos(ϕ)*sin(χ)
+#                    cos(ϕ)*cos(χ)-cos(θ)*sin(ϕ)*sin(χ)
+#                    sin(θ)*sin(χ)]
+#
+#     w = [sqrt(2/3) * 200 * m(54.735 * pi / 180, 45 * pi/180, x) for x=0:0.15:(2*pi+0.1)]
+#     base_vec =  [1 1 1] / sqrt(3)
+#
+#     for i=-5:5
+#         tt = [w[x] + vec(base_vec) + 50 * i for x=1:length(w)]
+#         x = map(x->tt[x][1], collect(1:length(w)))
+#         y = map(x->tt[x][2], collect(1:length(w)))
+#         z = map(x->tt[x][3], collect(1:length(w)))
+#         plot3D(x, y, z, color="blue")
+#     end
+#
+#     tt = [w[x] + vec(base_vec) + 50 * -5 for x=1:length(w)]
+#     x_start = map(x->tt[x][1], collect(1:length(w)))[1:5:end]
+#     y_start = map(x->tt[x][2], collect(1:length(w)))[1:5:end]
+#     z_start = map(x->tt[x][3], collect(1:length(w)))[1:5:end]
+#
+#
+#     tt = [w[x] + vec(base_vec) + 50 * 5 for x=1:length(w)]
+#     x_end = map(x->tt[x][1], collect(1:length(w)))[1:5:end]
+#     y_end = map(x->tt[x][2], collect(1:length(w)))[1:5:end]
+#     z_end = map(x->tt[x][3], collect(1:length(w)))[1:5:end]
+#
+#     for i=1:length(x_start)
+#         x = [x_start[i], x_end[i]]
+#         y = [y_start[i], y_end[i]]
+#         z = [z_start[i], z_end[i]]
+#         plot3D(x, y, z, color="blue")
+#     end
+#
+#
+#     info("Calculation finished")
+#     #PyPlot.plot(ee, ss)
+#     #=
+#     plot3D(eig_vals[:, 1], eig_vals[:, 2], eig_vals[:, 3], color="red")
+#     PyPlot.title("Stress path and von Mises yield surface")
+#     PyPlot.xlabel("Eig Stress 1")
+#     PyPlot.ylabel("Eig Stress 2")
+#     PyPlot.zlabel("Eig Stress 3")
+#     PyPlot.grid()
+#     PyPlot.show()
+#     =#
+# end
+
+#function test_von_mises_planestress_basic()
 
     steps = 1000
-    strain_max = 0.003
-    num_cycles = 3
-    E = 200.0e3
-    nu =  0.3
+    strain_max = 0.004
+    num_cycles = 1.
+    E = 200000.
+    nu = 0.3
     ν = 0.3
-    C = stiffnessTensor(E, ν)
-
-    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
-    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))
-
-    strain_last = zeros(Float64, (6))
-    strain_p = zeros(Float64, (6))
-    stress = zeros(Float64, (6, 1))
-    stress_y =  200.0
-    ss = Float64[]
-    ee = Float64[]
-
-    eig_stress = zeros(Float64, (3, 3))
-    eig_vals = zeros(Float64, (steps, 3))
-
-    function fill_tensor(a, b)
-        a[1, 1] = b[1]
-        a[2, 2] = b[2]
-        a[3, 3] = b[3]
-
-        a[1, 2] = b[6]
-        a[1, 3] = b[5]
-        a[2, 3] = b[4]
-
-        a[2, 1] = b[6]
-        a[3, 1] = b[5]
-        a[3, 2] = b[4]
-    end
-
-    mat = State(C, stress_y, zeros(Float64, 6), zeros(Float64, 6))
-
-    info("Starting calculation")
-    tic()
-    #=
-    for i=1:steps
-        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.stress)
-        eig_vals[i, :] = sort(eigvals(eig_stress))
-    end
-    =#
-    stress = zeros(Float64, 6)
-    strain = zeros(Float64, 6)
-    for i=1:steps
-        strain_new = reshape(strain_tot[i, :, :], (6, 1))
-        dstrain = strain_new - strain
-        calculate_stress!(dstrain, stress, C, stress_y, Val{:vonMises})
-        strain = vec(strain_new)
-        push!(ss, stress[1])
-        push!(ee, strain[1])
-        fill_tensor(eig_stress, stress)
-        eig_vals[i, :] = sort(eigvals(eig_stress))
-    end
-
-    toc()
-    # ================ Plotting =================== #
-    n(θ, ϕ) = [sin(θ)*cos(ϕ)
-               sin(θ)*sin(ϕ)
-               cos(θ)]
-    m(θ, ϕ, χ) = [-sin(ϕ)*cos(χ)-cos(θ)*cos(ϕ)*sin(χ)
-                   cos(ϕ)*cos(χ)-cos(θ)*sin(ϕ)*sin(χ)
-                   sin(θ)*sin(χ)]
-
-    w = [sqrt(2/3) * 200 * m(54.735 * pi / 180, 45 * pi/180, x) for x=0:0.15:(2*pi+0.1)]
-    base_vec =  [1 1 1] / sqrt(3)
-
-    for i=-5:5
-        tt = [w[x] + vec(base_vec) + 50 * i for x=1:length(w)]
-        x = map(x->tt[x][1], collect(1:length(w)))
-        y = map(x->tt[x][2], collect(1:length(w)))
-        z = map(x->tt[x][3], collect(1:length(w)))
-        plot3D(x, y, z, color="blue")
-    end
-
-    tt = [w[x] + vec(base_vec) + 50 * -5 for x=1:length(w)]
-    x_start = map(x->tt[x][1], collect(1:length(w)))[1:5:end]
-    y_start = map(x->tt[x][2], collect(1:length(w)))[1:5:end]
-    z_start = map(x->tt[x][3], collect(1:length(w)))[1:5:end]
-
-
-    tt = [w[x] + vec(base_vec) + 50 * 5 for x=1:length(w)]
-    x_end = map(x->tt[x][1], collect(1:length(w)))[1:5:end]
-    y_end = map(x->tt[x][2], collect(1:length(w)))[1:5:end]
-    z_end = map(x->tt[x][3], collect(1:length(w)))[1:5:end]
-
-    for i=1:length(x_start)
-        x = [x_start[i], x_end[i]]
-        y = [y_start[i], y_end[i]]
-        z = [z_start[i], z_end[i]]
-        plot3D(x, y, z, color="blue")
-    end
-
-
-    info("Calculation finished")
-    #PyPlot.plot(ee, ss)
-    #=
-    plot3D(eig_vals[:, 1], eig_vals[:, 2], eig_vals[:, 3], color="red")
-    PyPlot.title("Stress path and von Mises yield surface")
-    PyPlot.xlabel("Eig Stress 1")
-    PyPlot.ylabel("Eig Stress 2")
-    PyPlot.zlabel("Eig Stress 3")
-    PyPlot.grid()
-    PyPlot.show()
-    =#
-end
-
-function test_von_mises_planestress_basic()
-
-    steps = 1000
-    strain_max = 0.003
-    num_cycles = 5
-    E = 200.0e3
-    nu =  0.3
-    ν = 0.3
-    C = stiffnessTensorPlaneStress(E, ν)
+    C = E/((1+nu)*(1-2*nu)) .* [
+        1-nu   nu 0
+        nu   1-nu 0
+        0    0    (1-2*nu)/2]
 
     strain_tot = zeros(Float64, (steps, 3))
 
@@ -151,7 +155,7 @@ function test_von_mises_planestress_basic()
     strain_last = zeros(Float64, (3))
     strain_p = zeros(Float64, (3))
     stress = zeros(Float64, (3, 1))
-    stress_y =  200.0
+    stress_y =  400
     ss = Float64[]
     ee = Float64[]
 
@@ -161,35 +165,18 @@ function test_von_mises_planestress_basic()
 
     eig_stress = zeros(Float64, (3, 3))
     eig_vals = zeros(Float64, (steps, 3))
-    #mat = State(C, stress_y, zeros(Float64, 6), zeros(Float64, 6))
 
     info("Starting calculation")
     tic()
-    #=
-    for i=1:steps
-        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.stress)
-        eig_vals[i, :] = sort(eigvals(eig_stress))
-    end
-    =#
     stress = zeros(Float64, 3)
     strain = zeros(Float64, 3)
+    params = Dict("yield_stress" => stress_y)
+    Dtan = C
     for i=1:steps
         strain_new = reshape(strain_tot[i, :, :], (3, 1))
         dstrain = strain_new - strain
-        stress_inc, lambda = calculate_stress(dstrain,
-                                              stress,
-                                              C,
-                                              stress_y,
-                                              Val{:vonMises},
-                                              Val{:PlaneStressElasticPlasticProblem})
-        stress += stress_inc
+        JuliaFEM.plastic_von_mises!(stress, dstrain, C, params, Dtan)
         strain = vec(strain_new)
         s1, s2, t12  = stress
         se1 = (s1 + s2)/2 + sqrt(((s1 - s2)/2)^2 + t12^2)
@@ -197,14 +184,13 @@ function test_von_mises_planestress_basic()
         push!(ss, se1)
         push!(ee, se2)
     end
-
     toc()
 
     function vm_upper(a, c)
         vals = f(a[1], a[2], c)
         vm(vals[1], vals[2], 200)
     end
-    vm(a,b) = sqrt(a^2 - a*b + b^2) - 200
+    vm(a,b) = sqrt(a^2 - a*b + b^2) - stress_y
     f(m,c) = [600*cos(c) 600*sin(c)].*m
     x_vals = []
     max_iter = 100
@@ -229,15 +215,12 @@ function test_von_mises_planestress_basic()
         push!(x_vals, s11)
         push!(y_vals, s22)
     end
-    #=
-    PyPlot.plot(x_vals, y_vals)
-    PyPlot.plot(ee, ss)
-    PyPlot.grid()
-    PyPlot.show()
-    =#
-end
+
+    #plot(x_vals, y_vals)
+    plot(ee, ss)
+    show()
+# end
 
 # test_von_mises_3D_basic()
 
 #test_von_mises_planestress_basic()
-