diff --git a/src/JuliaFEM.jl b/src/JuliaFEM.jl
index e558b4a..a55fc49 100644
--- a/src/JuliaFEM.jl
+++ b/src/JuliaFEM.jl
@@ -19,7 +19,7 @@ include("api.jl")
 end
 
 module Preprocess
-include("abaqus_reader.jl") 
+include("abaqus_reader.jl")
 include("aster_reader.jl")
 end
 
@@ -32,6 +32,9 @@ module Test
 include("test.jl")
 end
 
+module MaterialModels
+include("vonmises.jl")
+end
 
 module Interfaces
 include("interfaces.jl")
diff --git a/src/vonMises.jl b/src/vonMises.jl
index c3bfd88..405b8ea 100644
--- a/src/vonMises.jl
+++ b/src/vonMises.jl
@@ -1,13 +1,17 @@
 # imports
+
+
+
 using ForwardDiff
 using NLsolve
-using PyPlot
+
 
 """
 Create a isotropic Hooke material matrix C
 
-More information: # http://www.efunda.com/formulae/solid_mechanics/mat_mechanics/hooke_isotropic.cfm
-
+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
@@ -19,6 +23,22 @@ 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)
+            end
+        end
+    end
+end
+=#
+#function outer_product(a, b)
+#
+#end
 function hookeStiffnessTensor(E, ν)
     a = 1 - ν
     b = 1 - 2*ν
@@ -33,10 +53,12 @@ function hookeStiffnessTensor(E, ν)
 end
 
 # Pick material values
-E = 200.0e3
-ν =  0.3
-C = hookeStiffnessTensor(E, ν)
-
+#E = 200.0e3
+#mu =  0.3
+#C = hookeStiffnessTensor(E, ν)
+#lambda =
+#nu =
+#C = λ * I ⊗ I + 2 * μ * II
 
 type State
     C :: Array{Float64, 2}
@@ -121,11 +143,11 @@ Returns
 -------
     Array{Float64, 7}, return values for solver
 """
-function G(params, dϵ, C, k, σ_begin)
+function vonMisesRoot(params, dϵ, C, σ_y, σ_begin)
 
     # Creating wrapper for gradient
-    yield(pars) = vonMisesYield(pars, k)
-    dfdσ = ForwardDiff.gradient(yield)
+    yield_wrap(pars) = vonMisesYield(pars, σ_y)
+    dfdσ = ForwardDiff.gradient(yield_wrap)
 
     # Stress rate
     dσ = params[1:6]
@@ -137,14 +159,14 @@ function G(params, dϵ, C, k, σ_begin)
 
     # Calculating equations
     function_1 = dσ - C * (dϵ - dϵp[1:6])
-    function_2 = yield(σ_tot)
+    function_2 = yield_wrap(σ_tot)
     [vec(function_1); function_2]
 end
 
 
 
 """
-Function which calculates the stress. Also handles if any yielding happens
+Stress for ideal plastic von Mises material model
 
 Parameters
 ----------
@@ -164,110 +186,26 @@ Returns
     Tuple
     Plastic strain rate dϵᵖ and new stress vector σ
 """
-function calculate_stress!(dϵ, mat::State)
+function calculate_stress!(dϵ, mat::State, ::Type{Val{:vonMises}})
     σ = mat.σ
     C = mat.C
     σ_y = mat.σ_y
     # Test stress
     σ_tria = σ + C * dϵ
 
-    # Calculating yield
+    # Calculating and checking for yield
     yield = vonMisesYield(σ_tria, σ_y)
-
     if yield > 0
         # Yielding happened
         # Creating functions for newton: xₙ₊₁ = xₙ - df⁻¹ * f and initial values
         initial_guess = [vec(σ_tria - σ); 0.1]
-        f(σ_)  = G(σ_, dϵ, C, σ_y, σ)
+        f(σ_)  = vonMisesRoot(σ_, dϵ, C, σ_y, σ)
         df     = ForwardDiff.jacobian(f)
 
         # Calculating root
         result = nlsolve(not_in_place(f, df), initial_guess).zero
-
         mat.σ  += result[1:6]
     else
         mat.σ = vec(σ_tria)
     end
 end
-
-steps = 1000
-strain_max = 0.003
-num_cycles = 3
-
-ϵ_tot = zeros(Float64, (steps, 6))
-ϵ_tot2 = zeros(Float64, (steps, 6))
-ϵ_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))
-
-ϵ_last = zeros(Float64, (6))
-ϵᵖ = zeros(Float64, (6))
-σ = zeros(Float64, (6, 1))
-σ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, σy, zeros(Float64, 6), zeros(Float64, 6))
-
-info("Starting calculation")
-for i=1:steps
-    ϵ_new = reshape(ϵ_tot[i, :, :], (6, 1))
-    dϵ = ϵ_new - mat.ϵ
-    calculate_stress!(dϵ, mat)
-    mat.ϵ += vec(dϵ)
-    push!(ss, mat.σ[1])
-    push!(ee, mat.ϵ[1])
-
-    fill_tensor(eig_stress, mat.σ)
-    eig_vals[i, :] = sort(eigvals(eig_stress))
-end
-
-# ================ 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.2:(2*pi+0.1)]
-base_vec =  [1 1 1] / sqrt(3)
-
-for i=-7:7
-    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
-
-
-#n(54.735 * pi / 180, 45 * pi/180)
-info("Calculation finished")
-#PyPlot.plot(ee, ss)
-plot3D(eig_vals[:, 1], eig_vals[:, 2], eig_vals[:, 3], color="red")
-PyPlot.title("Stress-Strain curve")
-PyPlot.xlabel("Strain")
-PyPlot.ylabel("Stress")
-PyPlot.grid()
-# PyPlot.plot(ee, ss)
-PyPlot.show()
diff --git a/test/test_von_mises_material.jl b/test/test_von_mises_material.jl
new file mode 100644
index 0000000..5829a81
--- /dev/null
+++ b/test/test_von_mises_material.jl
@@ -0,0 +1,123 @@
+module VonMisesTests
+
+using PyPlot
+
+macro R_str(s)
+    s
+end
+
+
+using JuliaFEM.MaterialModels: hookeStiffnessTensor, calculate_stress!, State
+
+function test_von_mises_basic()
+
+    steps = 1000
+    strain_max = 0.003
+    num_cycles = 3
+    E = 200.0e3
+    nu =  0.3
+    ν = 0.3
+    C = hookeStiffnessTensor(E, ν)
+
+    ϵ_tot = zeros(Float64, (steps, 6))
+    ϵ_tot2 = zeros(Float64, (steps, 6))
+    ϵ_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))
+
+    ϵ_last = zeros(Float64, (6))
+    ϵᵖ = zeros(Float64, (6))
+    σ = zeros(Float64, (6, 1))
+    σ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, σ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])
+
+        fill_tensor(eig_stress, mat.σ)
+        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
+
+test_von_mises_basic()
+
+end