diff --git a/src/elasticity.jl b/src/elasticity.jl
index 62890ba..8275fb0 100644
--- a/src/elasticity.jl
+++ b/src/elasticity.jl
@@ -1,13 +1,11 @@
 # This file is a part of JuliaFEM.
 # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
 
-include("vonmises.jl")
+include("elasticplastic.jl")
 
 # Elasticity problems
-
-abstract ElasticityProblem     <: AbstractProblem
-
-abstract PlaneStressElasticityProblem     <: ElasticityProblem
+abstract ElasticityProblem            <: AbstractProblem
+abstract PlaneStressElasticityProblem <: ElasticityProblem
 
 function get_unknown_field_name{P<:ElasticityProblem}(::Type{P})
     return "displacement"
@@ -17,14 +15,6 @@ 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
-
 # 3D Elasticity problems
 function ElasticityProblem(dim::Int=3, elements=[])
     return Problem{ElasticityProblem}("elasticity problem", dim, elements)
@@ -67,8 +57,6 @@ https://en.wikipedia.org/wiki/Hooke's_law
 
 """
 function get_residual_vector{P<:ElasticityProblem}(problem::Problem{P}, element::Element, ip::IntegrationPoint, time::Number; variation=nothing)
-#function get_residual_vector{P<:ElasticPlasticProblem}(problem::Problem{P}, element::Element, ip::IntegrationPoint, time::Number; variation=nothing)
-
     r = zeros(Float64, problem.dim, length(element))
 
     J = get_jacobian(element, ip, time)
diff --git a/src/elasticplastic.jl b/src/elasticplastic.jl
new file mode 100644
index 0000000..a919df0
--- /dev/null
+++ b/src/elasticplastic.jl
@@ -0,0 +1,247 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
+
+include("vonmises.jl")
+
+# Elasticity problems
+abstract ElasticPlasticProblem            <: AbstractProblem
+abstract PlaneStressElasticPlasticProblem <: ElasticPlasticProblem
+
+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
+
+# 3D Elasticity problems
+function ElasticPlasticProblem(dim::Int=3, elements=[])
+    return Problem{ElasticPlasticProblem}("elasticplastic problem", dim, elements)
+end
+
+# 2D Plane stress elasticity problems
+function PlaneStressElasticPlasticProblem(dim::Int=2, elements=[])
+    return Problem{PlaneStressElasticPlasticProblem}("plane stress elasticplastic problem", dim, elements)
+end
+
+
+function get_residual_vector{P<:PlaneStressElasticPlasticProblem}(problem::Problem{P}, element::Element, ip::IntegrationPoint, time::Number; variation=nothing)
+    r = zeros(Float64, problem.dim, length(element))
+
+    J = get_jacobian(element, ip, time)
+
+    # internal forces
+    if haskey(element, "youngs modulus") && haskey(element, "poissons ratio")
+        if !haskey(element, "integration points")
+            if P == PlaneStressElasticPlasticProblem
+                last_stress = zeros(2,2)
+                last_strain = zeros(2,2)
+            else
+                last_stress = zeros(3,3)
+                last_strain = zeros(3,3)
+            end
+        else
+            for each_ip in element("integration points", time)
+                if isapprox(each_ip.xi, ip.xi)
+                    last_stress = ip("stress", time)
+                    last_strain = ip("stress", time)
+                    break
+                end
+            end
+        end
+        u = element("displacement", time, variation)
+        grad = element(ip, time, Val{:grad})
+        gradu = grad*u
+
+        # deformation gradient
+        F = I + gradu
+        E = 1/2*(F'*F - I)
+
+        # material
+        young = element("youngs modulus", ip, time)
+        poisson = element("poissons ratio", ip, time)
+        stress_y = element("yield stress", time).data
+        mu = young/(2*(1+poisson))
+        lambda = young*poisson/((1+poisson)*(1-2*poisson))
+        if P == PlaneStressElasticPlasticProblem
+            lambda = 2*lambda*mu/(lambda + 2*mu)  # <- correction for 2d problems
+        end
+        dstrain = E - last_strain
+        material_model = element("material model", time)
+        s = last_stress
+        de = copy(ForwardDiff.get_value(dstrain))
+
+        if P == PlaneStressElasticPlasticProblem
+            C = stiffnessTensorPlaneStress(young, poisson)
+            s_v = [s[1,1], s[2,2], s[1,2]]
+            de_ = [de[1,1], de[2,2], de[1,2]]
+            problem_stress_type = :PlaneStressElasticPlasticProblem
+        else
+            C = stiffnessTensor(young, poisson)
+            s_v = [s[1,1], s[2,2], s[3,3], s[2,3], s[1,3], s[1,2]]
+            de_ = [de[1,1], de[2,2], de[3,3], de[2,3], de[1,3], de[1,2]]
+            problem_stress_type = :ElasticPlasticProblem
+        end
+        dep = zeros(3)
+        stress_inc, dep = calculate_stress(de_,
+                                           s_v,
+                                           C,
+                                           stress_y,
+                                           Val{:vonMises},
+                                           Val{problem_stress_type})
+        nd = [dep[1] dep[3];
+              dep[3] dep[2]]
+    #    a = dstrain - nd
+    #    b = C * a
+        info("%% ", dep)
+        dif = dstrain - nd
+        mm = [dif[1,1], dif[2,2], dif[1,2]]
+        s_v += C * mm
+        info("--: ", ForwardDiff.get_value(s_v))
+        # stress
+        if P == PlaneStressElasticPlasticProblem
+            S = [s_v[1] s_v[3];
+                 s_v[3] s_v[2]]
+        else
+        S = [s_v[1] s_v[6] s_v[5];
+             s_v[6] s_v[2] s_v[4];
+             s_v[5] s_v[4] s_v[3]]
+        end
+        r += F*S*grad*det(J)
+    end
+
+    # external forces - volume load
+    if haskey(element, "displacement load")
+        basis = element(ip, time)
+        b = element("displacement load", ip, time)
+        r -= b*basis*det(J)
+    end
+
+    # external forces - surface traction force
+    if haskey(element, "displacement traction force")
+        basis = element(ip, time)
+        T = element("displacement traction force", ip, time)
+        JT = transpose(J)
+        s = size(JT, 2) == 1 ? JT : cross(JT[:,1], JT[:,2])
+        r -= T*basis*norm(s)
+    end
+
+    return vec(r)
+end
+
+
+
+#=
+function get_residual_vector{P<:ElasticPlasticProblem}(problem::Problem{P}, element::Element, ip::IntegrationPoint, time::Number; variation=nothing)
+    r = zeros(Float64, problem.dim, length(element))
+
+    J = get_jacobian(element, ip, time)
+
+    info("_____________________")
+    # internal forces
+    if haskey(element, "youngs modulus") && haskey(element, "poissons ratio")
+
+        if !haskey(element, "integration points")
+            if P == PlaneStressElasticPlasticProblem
+                last_stress = zeros(2,2)
+                last_strain = zeros(2,2)
+            else
+                last_stress = zeros(3,3)
+                last_strain = zeros(3,3)
+            end
+        else
+            for each_ip in element("integration points", time)
+                if isapprox(each_ip.xi, ip.xi)
+                    last_stress = ip("stress", time)
+                    last_strain = ip("stress", time)
+                    break
+                end
+            end
+        end
+        u = element("displacement", time, variation)
+        grad = element(ip, time, Val{:grad})
+        gradu = grad*u
+
+        # deformation gradient
+        F = I + gradu
+
+        # material
+        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
+
+        # strain
+        E = 1/2*(F'*F - I)
+        #E = 1/2*(gradu + gradu') # finite strain (total)
+
+        young = element("youngs modulus", ip, time)
+        poisson = element("poissons ratio", ip, time)
+        stress_y = element("yield stress", time).data
+        dstrain = E - last_strain
+        material_model = element("material model", time)
+        s = last_stress
+        de = ForwardDiff.get_value(dstrain)
+
+        if P == PlaneStressElasticPlasticProblem
+            C = stiffnessTensorPlaneStress(young, poisson)
+            s_v = [s[1,1], s[2,2], s[1,2]]
+            de_ = [de[1,1], de[2,2], de[1,2]]
+            problem_stress_type = :PlaneStressElasticPlasticProblem
+        else
+            C = stiffnessTensor(young, poisson)
+            s_v = [s[1,1], s[2,2], s[3,3], s[2,3], s[1,3], s[1,2]]
+            de_ = [de[1,1], de[2,2], de[3,3], de[2,3], de[1,3], de[1,2]]
+            problem_stress_type = :ElasticPlasticProblem
+        end
+
+        stress_inc, lambda = plastic_multiplier = calculate_stress(de_,
+                                                           s_v,
+                                                           C,
+                                                           stress_y,
+                                                           Val{:vonMises},
+                                                           Val{problem_stress_type})
+
+        # dep = lambda * dfds(s)
+        # upate_material_parameters!(...)
+        s_new = s_v + stress_inc
+        #S = [s_v[1] s_v[6] s_v[5];
+    #         s_v[6] s_v[2] s_v[4];
+#             s_v[5] s_v[4] s_v[3]]
+        S = [s_new[1] s_new[3];
+             s_new[3] s_new[2]]
+        # S = C * (E - dep)
+
+
+        info("Stress: ", vec(ForwardDiff.get_value(S)))
+        # stress
+        #S = lambda*trace(E)*I + 2*mu*E
+
+        r += F*S*grad*det(J)
+
+    end
+
+
+    # external forces - volume load
+    if haskey(element, "displacement load")
+        basis = element(ip, time)
+        b = element("displacement load", ip, time)
+        r -= b*basis*det(J)
+    end
+
+    # external forces - surface traction force
+    if haskey(element, "displacement traction force")
+        basis = element(ip, time)
+        T = element("displacement traction force", ip, time)
+        JT = transpose(J)
+        s = size(JT, 2) == 1 ? JT : cross(JT[:,1], JT[:,2])
+        r -= T*basis*norm(s)
+    end
+
+    return vec(r)
+end
+=# #fff
diff --git a/src/vonmises.jl b/src/vonmises.jl
index f032713..653bd67 100644
--- a/src/vonmises.jl
+++ b/src/vonmises.jl
@@ -1,5 +1,4 @@
 using ForwardDiff
-# using NLsolve
 
 """
 Create a isotropic Hooke material matrix C
@@ -31,6 +30,17 @@ function stiffnessTensor(E, ν)
                    0 0 0 0 0 b].*multiplier
 end
 
+# Creating functions for newton: xₙ₊₁ = xₙ - df⁻¹ * f and initial values
+function find_root!(f, df, x; max_iter=50, norm_acc=1e-10)
+    converged = false
+    for i=1:max_iter
+        dx = df(x) \ -f(x)
+        x += dx
+        norm(dx) < norm_acc && (converged = true; break)
+    end
+    converged || error("no convergence!")
+    x
+end
 
 type State
     C :: Array{Float64, 2}
@@ -176,7 +186,9 @@ function calculate_stress!(dstrain, mat::State, ::Type{Val{:vonMises}})
     end
 end
 
-function calculate_stress!(dstrain, stress, C, stress_y, ::Type{Val{:vonMises}})
+function calculate_stress(dstrain, stress, C, stress_y,
+     ::Type{Val{:vonMises}},
+     ::Type{Val{:ElasticPlasticProblem}})
     # Test stress
     stress_tria = stress + C * dstrain
 
@@ -184,7 +196,7 @@ function calculate_stress!(dstrain, stress, C, stress_y, ::Type{Val{:vonMises}})
     yield = vonMisesYield(stress_tria, stress_y)
     if isless(yield, 0.0)
         # stress[i] = stress_tria[i]
-        return 0.0 
+        return 0.0
     else
         # Yielding happened
         # Creating functions for newton: xₙ₊₁ = xₙ - df⁻¹ * f and initial values
@@ -193,7 +205,7 @@ function calculate_stress!(dstrain, stress, C, stress_y, ::Type{Val{:vonMises}})
         df     = ForwardDiff.jacobian(f)
 
         # Calculating root
-        # result = nlsolve(not_in_place(f, df), initial_guess).zero 
+        # result = nlsolve(not_in_place(f, df), initial_guess).zero
         max_iter = 10
         converged = false
         for i=1:5
@@ -208,3 +220,83 @@ function calculate_stress!(dstrain, stress, C, stress_y, ::Type{Val{:vonMises}})
     end
 end
 
+##################################################################################
+# ----- AFTER THIS POINT: VON MISES : PLANE STRESS IMPLEMENTATION          ----- #
+##################################################################################
+
+"""
+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)
+    s1, s2, t12  = stress
+    # Calculating principal stresses
+    # 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
+
+# https://andriandriyana.files.wordpress.com/2008/03/yield_criteria.pdf
+function vonMisesYieldPlaneStress(stress, stress_y)
+    stress_eq_plane_stress(stress) - stress_y
+end
+
+function vonMisesRootPlaneStress(params, dstrain, C, stress_y, stress_base)
+
+    # Creating wrapper for gradient
+    vm_wrap(stress_) = vonMisesYieldPlaneStress(stress_, stress_y)
+    dfds = ForwardDiff.gradient(vm_wrap)
+
+    # Stress rate and total strain
+    dstress = params[1:3]
+    stress_tot = vec(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_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}})
+    # Test stress
+    dstress = C * dstrain
+    stress_tria = stress + dstress
+
+    # Calculating and checking for yield
+    yield = vonMisesYieldPlaneStress(stress_tria, stress_y)
+    if isless(yield, 0.0)
+        return dstress, zeros(3)
+    else
+        info("yielded")
+        # Yielding happened
+        # 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)
+        # Calculating root
+        results = find_root!(f, df, x)
+        stress_tot = stress + results[1:3]
+        plastic_multiplier = results[end]
+        vm_wrap(stress_) = vonMisesYieldPlaneStress(stress_, stress_y)
+        dfds = ForwardDiff.gradient(vm_wrap)
+        dep = plastic_multiplier * dfds(stress_tot)
+        info("II ", stress_tot)
+        return results[1:3], dep
+    end
+end
diff --git a/test/test_directsolver_with_vonmises.jl b/test/test_directsolver_with_vonmises.jl
new file mode 100644
index 0000000..a044cc8
--- /dev/null
+++ b/test/test_directsolver_with_vonmises.jl
@@ -0,0 +1,230 @@
+
+
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
+
+module DirectSolverVonMisesTests
+
+using JuliaFEM.Test
+using JuliaFEM.Core: Seg2, Quad4
+using JuliaFEM.Core: PlaneStressElasticityProblem, DirichletProblem
+using JuliaFEM.Core: PlaneStressElasticPlasticProblem
+using JuliaFEM.Core: DirectSolver
+
+function test_solver_multiple_dirichlet_bc()
+
+    N = Vector[[0.0, 0.0], [1.0, 0.0], [0.0, 1.0], [1.0, 1.0]]
+
+    e1 = Quad4([1, 2, 4, 3])
+    e1["geometry"] = Vector[N[1], N[2], N[4], N[3]]
+    e1["youngs modulus"] = 900.0
+    e1["poissons ratio"] = 0.25
+    e1["yield stress"] = 100.0
+    e1["material model"] = :vonMises
+    b1 = Seg2([3, 4])
+    b1["geometry"] = Vector[N[3], N[4]]
+    b1["displacement traction force"] = (
+        0.0 => Vector[[0.0,    0.0], [0.0,    0.0]],
+        1.0 => Vector[[0.0, -100.0], [0.0, -100.0]])
+
+    #problem = PlaneStressElasticityProblem()
+    problem = PlaneStressElasticPlasticProblem()
+    push!(problem, e1)
+    push!(problem, b1)
+
+    # boundary elements for dirichlet dx=0
+    dx = Seg2([1, 3])
+    dx["geometry"] = Vector[N[1], N[3]]
+    dx["displacement 1"] = 0.0
+
+    # boundary elements for dirichlet dy=0
+    dy = Seg2([1, 2])
+    dy["geometry"] = Vector[N[1], N[2]]
+    dy["displacement 2"] = 0.0
+
+    problem2 = DirichletProblem("displacement", 2)
+    push!(problem2, dx)
+
+    problem3 = DirichletProblem("displacement", 2)
+    push!(problem3, dy)
+
+    solver = DirectSolver()
+    #solver.dump_matrices = true
+    solver.name = "test_solver_multiple_dirichlet_bc"
+    push!(solver, problem)
+    push!(solver, problem2)
+    push!(solver, problem3)
+
+    # launch solver
+    #norm = solver(0.0)
+    norm = solver(1.0)
+    disp = e1("displacement", [1.0, 1.0], 1.0)
+    info("displacement at tip: $disp")
+    #@test isapprox(disp, [3.17431158889468E-02, -1.38591518927826E-01])
+
+end
+test_solver_multiple_dirichlet_bc()
+
+#=
+function test_direct_cholesky_with_non_homogeneous_dirichlet_conditions()
+
+    N = Vector[[0.0, 0.0], [1.0, 0.0], [0.0, 1.0], [1.0, 1.0]]
+
+    e1 = Quad4([1, 2, 4, 3])
+    e1["geometry"] = Vector[N[1], N[2], N[4], N[3]]
+    e1["youngs modulus"] = 900.0
+    e1["poissons ratio"] = 0.25
+
+    problem = PlaneStressElasticityProblem()
+    push!(problem, e1)
+
+    # left boundary: dx=-0.1, dy=0.1
+    bc1 = Seg2([1, 3])
+    bc1["geometry"] = Vector[N[1], N[3]]
+    bc1["displacement 1"] = -0.1
+    bc1["displacement 2"] =  0.1
+
+    # right boundary: dx=0.2, dy=-0.2
+    bc2 = Seg2([2, 4])
+    bc2["geometry"] = Vector[N[2], N[4]]
+    bc2["displacement 1"] =  0.2
+    bc2["displacement 2"] = -0.2
+
+    boundary = DirichletProblem("displacement", 2)
+    push!(boundary, bc1)
+    push!(boundary, bc2)
+
+    solver = DirectSolver("test_direct_cholesky_with_non_homogeneous_dirichlet_boundary_conditions")
+    push!(solver, problem)
+    push!(solver, boundary)
+
+    # launch solver
+    solver.method = :UMFPACK
+    solver.dump_matrices = true
+    solver.max_iterations = 1
+    iters, status = solver(0.0)
+#   FIXME: solver gives no convergence warning when all dofs are fixed.
+    n1disp = e1("displacement", [-1.0, -1.0], 0.0)
+    n2disp = e1("displacement", [ 1.0, -1.0], 0.0)
+    n3disp = e1("displacement", [-1.0,  1.0], 0.0)
+    n4disp = e1("displacement", [ 1.0,  1.0], 0.0)
+    udisp = [n1disp n2disp n3disp n4disp]
+    info("nodal disp = ", udisp)
+    @test isapprox(n1disp, [-0.1,  0.1])
+    @test isapprox(n3disp, [-0.1,  0.1])
+    @test isapprox(n2disp, [ 0.2, -0.2])
+    @test isapprox(n4disp, [ 0.2, -0.2])
+    @test status == true
+end
+#test_direct_cholesky_with_non_homogeneous_dirichlet_conditions()
+
+function test_solver_no_convergence()
+
+    N = Vector[[0.0, 0.0], [1.0, 0.0], [0.0, 1.0], [1.0, 1.0]]
+
+    e1 = Quad4([1, 2, 4, 3])
+    e1["geometry"] = Vector[N[1], N[2], N[4], N[3]]
+    e1["youngs modulus"] = 900.0
+    e1["poissons ratio"] = 0.25
+    b1 = Seg2([3, 4])
+    b1["geometry"] = Vector[N[3], N[4]]
+    b1["displacement traction force"] = Vector[[100.0, 100.0], [100.0, 100.0]]
+
+    problem = PlaneStressElasticityProblem()
+    push!(problem, e1)
+    push!(problem, b1)
+
+    # boundary elements for dirichlet dx=0
+    dx = Seg2([1, 3])
+    dx["geometry"] = Vector[N[1], N[3]]
+    dx["displacement 1"] = 0.0
+
+    # boundary elements for dirichlet dy=0
+    dy = Seg2([1, 2])
+    dy["geometry"] = Vector[N[1], N[2]]
+    dy["displacement 2"] = 0.0
+
+    problem2 = DirichletProblem("displacement", 2)
+    push!(problem2, dx)
+
+    problem3 = DirichletProblem("displacement", 2)
+    push!(problem3, dy)
+
+    solver = DirectSolver()
+    solver.max_iterations = 1
+    push!(solver, problem)
+    push!(solver, problem2)
+    push!(solver, problem3)
+
+    # launch solver
+    iterations, status = solver(0.0)
+    @test status == false
+end
+
+
+function test_solver_multiple_bodies_multiple_dirichlet_bc()
+    N = Vector[
+        [0.0, 0.0], [1.0, 0.0],
+        [0.0, 1.0], [1.0, 1.0],
+        [0.0, 2.0], [1.0, 2.0]]
+
+    e1 = Quad4([1, 2, 4, 3])
+    e1["geometry"] = Vector[N[1], N[2], N[4], N[3]]
+    e2 = Quad4([3, 4, 6, 5])
+    e2["geometry"] = Vector[N[3], N[4], N[6], N[5]]
+    for el in [e1, e2]
+        el["youngs modulus"] = 900.0
+        el["poissons ratio"] = 0.25
+    end
+    b1 = Seg2([5, 6])
+    b1["geometry"] = Vector[N[5], N[6]]
+    b1["displacement traction force"] = Vector[[0.0, -100.0], [0.0, -100.0]]
+
+    body1 = PlaneStressElasticityProblem()
+    push!(body1, e1)
+
+    body2 = PlaneStressElasticityProblem()
+    push!(body2, e2)
+    push!(body2, b1)
+
+    # boundary elements for dirichlet dx=0
+    dx1 = Seg2([1, 3])
+    dx1["geometry"] = Vector[N[1], N[3]]
+    dx2 = Seg2([3, 5])
+    dx2["geometry"] = Vector[N[3], N[5]]
+    for dx in [dx1, dx2]
+        dx["displacement 1"] = 0.0
+    end
+
+    boundary1 = DirichletProblem("displacement", 2)
+    push!(boundary1, dx1)
+    push!(boundary1, dx2)
+
+    # boundary elements for dirichlet dy=0
+    dy1 = Seg2([1, 2])
+    dy1["geometry"] = Vector[N[1], N[2]]
+    dy1["displacement 2"] = 0.0
+
+    boundary2 = DirichletProblem("displacement", 2)
+    push!(boundary2, dy1)
+
+
+    solver = DirectSolver()
+    push!(solver, body1)
+    push!(solver, body2)
+    push!(solver, boundary1)
+    push!(solver, boundary2)
+
+    # launch solver
+    norm = solver(0.0)
+
+    disp = e2("displacement", [1.0, 1.0], 0.0)
+    info("displacement at tip: $disp")
+    # code aster verification, two_elements.comm
+    @test isapprox(disp, [3.17431158889468E-02, -2.77183037855653E-01])
+
+end
+
+#test_solver_multiple_bodies_multiple_dirichlet_bc()
+=#
+end
diff --git a/test/test_von_mises_material.jl b/test/test_von_mises_material.jl
index 9ece1a4..c696883 100644
--- a/test/test_von_mises_material.jl
+++ b/test/test_von_mises_material.jl
@@ -1,9 +1,12 @@
 module VonMisesTests
 
 using PyPlot
-using JuliaFEM.MaterialModels: stiffnessTensor, calculate_stress!, State
+using JuliaFEM.Test
+using JuliaFEM.MaterialModels: stiffnessTensor, calculate_stress, State
+using JuliaFEM.MaterialModels: stiffnessTensorPlaneStress
 
-function test_von_mises_basic()
+
+function test_von_mises_3D_basic()
 
     steps = 1000
     strain_max = 0.003
@@ -51,7 +54,7 @@ function test_von_mises_basic()
 
     info("Starting calculation")
     tic()
-    #= 
+    #=
     for i=1:steps
         strain_new = reshape(strain_tot[i, :, :], (6, 1))
         dstrain = strain_new - mat.strain
@@ -76,7 +79,7 @@ function test_von_mises_basic()
         fill_tensor(eig_stress, stress)
         eig_vals[i, :] = sort(eigvals(eig_stress))
     end
-    
+
     toc()
     # ================ Plotting =================== #
     n(θ, ϕ) = [sin(θ)*cos(ϕ)
@@ -127,6 +130,112 @@ function test_von_mises_basic()
     PyPlot.show()
 end
 
-test_von_mises_basic()
+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, ν)
+
+    strain_tot = zeros(Float64, (steps, 3))
+
+    # 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_last = zeros(Float64, (3))
+    strain_p = zeros(Float64, (3))
+    stress = zeros(Float64, (3, 1))
+    stress_y =  200.0
+    ss = Float64[]
+    ee = Float64[]
+
+
+    ss2 = Float64[]
+    ee2 = Float64[]
+
+    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)
+    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
+        strain = vec(strain_new)
+        s1, s2, t12  = stress
+        se1 = (s1 + s2)/2 + sqrt(((s1 - s2)/2)^2 + t12^2)
+        se2 = (s1 + s2)/2 - sqrt(((s1 - s2)/2)^2 + t12^2)
+        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
+    f(m,c) = [600*cos(c) 600*sin(c)].*m
+    x_vals = []
+    max_iter = 100
+    y_vals = []
+    for i=0:0.1:(2*pi+0.3)
+        wf(x) = f(x, i)
+        t = 0.01
+        step = 2
+        merkki = -1
+        s11, s22 = wf(t)
+        ii = 0
+        while (abs(vm(s11, s22)) > 1e-7) && ii < max_iter
+            val = vm(s11, s22)
+            if sign(val) != merkki
+                merkki *= -1
+                step *= -0.5
+            end
+            t += step
+            s11, s22 = wf(t)
+            ii += 1
+        end
+        push!(x_vals, s11)
+        push!(y_vals, s22)
+    end
+    PyPlot.plot(x_vals, y_vals)
+    PyPlot.plot(ee, ss)
+    PyPlot.grid()
+    PyPlot.show()
+end
+
+# test_von_mises_3D_basic()
+
+test_von_mises_planestress_basic()
 
 end