diff --git a/test/test_elasticity_solver.jl b/test/test_elasticity_solver.jl
index 142eaf3..167c9e3 100644
--- a/test/test_elasticity_solver.jl
+++ b/test/test_elasticity_solver.jl
@@ -280,3 +280,122 @@ facts("test that elimination of non-homogeneous dirichlet boundary conditions ra
     I, J, V = findnz(A)
     @fact_throws I, V = eliminate_boundary_conditions(dirichletbc, I, V)
 end
+
+module TestElasticitySolver
+
+using JuliaFEM.elasticity_solver: calc_local_matrices
+
+facts("test solve one element model") do
+    X = [0.0 0.0; 10.0 0.0; 10.0 1.0; 0.0 1.0]'
+    F = [0 0; 0 0; 0 -2; 0 0]'
+
+    # Material properties
+    E = 90
+    nu = 0.25
+    mu = E/(2*(1+nu))
+    la = E*nu/((1+nu)*(1-2*nu))
+    la = 2*la*mu/(la + 2*mu)
+
+    u = zeros(2, 4)
+    du = zeros(2, 4)
+    R = zeros(2, 4)
+    K = zeros(8, 8)
+
+    basis(xi) = [
+        (1-xi[1])*(1-xi[2])/4
+        (1+xi[1])*(1-xi[2])/4
+        (1+xi[1])*(1+xi[2])/4
+        (1-xi[1])*(1+xi[2])/4]
+
+    dbasis(xi) = [-(1-xi[2])/4.0    -(1-xi[1])/4.0
+                   (1-xi[2])/4.0    -(1+xi[1])/4.0
+                   (1+xi[2])/4.0     (1+xi[1])/4.0
+                  -(1+xi[2])/4.0     (1-xi[1])/4.0]
+
+    ipoints = 1/sqrt(3)*[-1 -1; 1 -1; 1 1; -1 1]
+    iweights = [1, 1, 1, 1]
+    free_dofs = [3, 4, 5, 6]
+
+    for i=1:10
+        calc_local_matrices!(X, u, R, K, basis, dbasis, la, mu, ipoints, iweights)
+        du[free_dofs] = K[free_dofs, free_dofs] \ -(R - F)[free_dofs]
+        u += du
+        if norm(du) < 1.0e-9
+            Logging.debug("Converged in $i iterations.")
+            break
+        end
+    end
+
+    # Tested against Elmer solution
+    Logging.debug("solution vector: \n $u")
+    @fact u[2, 3] --> roughly(-2.222244754401764)
+    norm1 = norm(u)
+    Logging.debug("norm of u: $(norm(u))")
+
+    # We rotate model a bit and make sure that L2 norm is same
+    phi = 30/180*pi
+    rmat = [
+        cos(phi) -sin(phi)
+        sin(phi)  cos(phi)]
+    X = rmat*X
+    F = rmat*F
+    u = zeros(2, 4)
+    for i=1:10
+        calc_local_matrices!(X, u, R, K, basis, dbasis, la, mu, ipoints, iweights)
+        du[free_dofs] = K[free_dofs, free_dofs] \ -(R - F)[free_dofs]
+        u += du
+        if norm(du) < 1.0e-9
+            Logging.debug("Converged in $i iterations.")
+            break
+        end
+    end
+    Logging.debug("solution vector: \n $u")
+    Logging.debug("norm of u: $(norm(u))")
+    @fact norm(u) --> roughly(norm1) 
+
+    # test two element model
+    X = [0.0 0.0; 5.0 0.0; 5.0 1.0; 0.0 1.0]'
+    u = zeros(2, 6)
+    du = zeros(2, 6)
+    R = zeros(2, 4)
+    K = zeros(8, 8)
+    ass1 = [9, 10, 1, 2, 5, 6, 11, 12]
+    ass2 = [1, 2, 3, 4, 7, 8, 5, 6]
+    free_dofs = collect(1:8)
+    F = [0 0; 0 0; 0 0; 0 -0.1; 0 0; 0 0]'
+
+    A = zeros(12, 12)
+    b = zeros(2, 6)
+    for i=1:1
+        Logging.debug("Iteration $i")
+        A[:,:] = 0.0
+        b[:] = 0.0
+        #Logging.debug("Assembling")
+        for ass in (ass1, ass2)
+            #Logging.debug("ass = $ass, u[ass] = $(u[ass])")
+            calc_local_matrices!(X, u[ass], R, K, basis, dbasis, la, mu, ipoints, iweights)
+            A[ass,ass] += K
+            b[ass] += R[:]
+        end
+        dump(round(A, 2))
+        println("K norm = $(norm(A[free_dofs, free_dofs]))")
+        du[free_dofs] = A[free_dofs, free_dofs] \ -(b - F)[free_dofs]
+        println("du = $du")
+        u += du
+        Logging.debug("Norm of du: $(norm(du))")
+        for ass in (ass1, ass2)
+            Logging.debug("Element displacement: $(reshape(u[ass], 2, 4))")
+        end
+        if norm(du) < 1.0e-9
+            Logging.debug("Converged in $i iterations.")
+            break
+        end
+    end
+    Logging.debug("solution vector: \n $u")
+    Logging.debug("norm of u: $(norm(u))")
+    @pending norm(u) --> :something
+end
+
+exitstatus()
+
+end