diff --git a/src/elasticity_solver.jl b/src/elasticity_solver.jl
index d8a5b11..f9a429b 100644
--- a/src/elasticity_solver.jl
+++ b/src/elasticity_solver.jl
@@ -2,7 +2,7 @@
 module elasticity_solver
 
 using Logging
-@Logging.configure(level=DEBUG)
+@Logging.configure(level=INFO)
 
 VERSION < v"0.4-" && using Docile
 
@@ -178,8 +178,6 @@ function assemble!(fe, eldofs_, I, V)
                 push!(eldofs, dim*(i-1)+d)
             end
         end
-        @debug("old eldofs", eldofs_)
-        @debug("new eldofs", eldofs)
     end
 
     for i in 1:n
@@ -283,40 +281,29 @@ function solve_elasticity_increment!(X, u, du, elmap, nodalloads,
         elmap = elmap''
     end
     nelnodes, nelements = size(elmap)
-    #@debug("elements in model: ", nelements)
-    #@debug("nodes / element in model: ", nelnodes)
+    dim, nnodes = size(u)
+    dofs = dim*nelnodes
 
-
-    dim = size(u)[1]
     Imat = Int64[]
     Jmat = Int64[]
     Vmat = Float64[]
     Ivec = Int64[]
     Vvec = Float64[]
-    # FIXME: different number of elements / node
-    #@debug("problem dimension: ", dim)
 
-    dofs = dim*nelnodes
+    # FIXME: different number of nodes/element
     R = zeros(dim, nelnodes)
     Kt = zeros(dofs, dofs)
-    #@debug("size of R: ", size(R))
-    #@debug("size of Kt: ", size(Kt))
 
     # this can be parallelized
     for i in 1:nelements
         eldofs = elmap[:,i]
-        #@debug("element dofs ", eldofs)
-        #@debug("Assembling element ", i)
-        #@debug("Local coords:\n", X[:, eldofs])
-        #@debug("Local u:\n", u[:, eldofs])
-        calc_local_matrices!(X[:, eldofs], u[:, eldofs], R, Kt, N, dNdξ, λ, μ, ipoints, iweights)
-        #@debug("Assemble Kt")
+        calc_local_matrices!(X[:, eldofs], u[:, eldofs], R, Kt, N, dNdξ,
+                             λ[eldofs], μ[eldofs], ipoints, iweights)
         assemble!(Kt, eldofs, Imat, Jmat, Vmat)
-        #@debug("Assemble R")
         assemble!(R, eldofs, Ivec, Vvec)
     end
 
-    # add additional neumann boundary conditions
+    # add additional neumann boundary conditions to force vector
     for (i, nodal_load) in enumerate(nodalloads)
         if nodal_load == 0
             continue
@@ -325,19 +312,19 @@ function solve_elasticity_increment!(X, u, du, elmap, nodalloads,
         push!(Vvec, -nodal_load)
     end
 
-    # Remove dirichlet boundary conditions
-    Imat, Jmat, Vmat = eliminate_boundary_conditions(dirichletbc, Imat, Jmat, Vmat)
-    Ivec, Vvec = eliminate_boundary_conditions(dirichletbc, Ivec, Vvec)
-  #R -= nodalloads
-
-    # Create sparse matrices
+    # Create sparse matrix and vector
     A = sparse(Imat, Jmat, Vmat)
     b = sparsevec(Ivec, Vvec)
 
-    # solution
+    # Remove dirichlet boundary conditions
     free_dofs = find(isnan(dirichletbc))
-    #du[free_dofs] = Kt[free_dofs, free_dofs] \ -reshape(R, 8)[free_dofs]
+    #Imat, Jmat, Vmat = eliminate_boundary_conditions(dirichletbc, Imat, Jmat, Vmat)
+    b = b[free_dofs]
+    A = A[free_dofs, free_dofs]
+
+    # solution
     du[free_dofs] = lufact(A) \ -full(b)
 end
 
+
 end
diff --git a/test/test_elasticity_solver.jl b/test/test_elasticity_solver.jl
index 76a1fc7..cc29d23 100644
--- a/test/test_elasticity_solver.jl
+++ b/test/test_elasticity_solver.jl
@@ -1,6 +1,6 @@
 using FactCheck
 using Logging
-@Logging.configure(level=DEBUG)
+@Logging.configure(level=INFO)
 
 
 using JuliaFEM.elasticity_solver: solve_elasticity_increment!
@@ -86,6 +86,57 @@ facts("test solve elasticity increment rot 30") do
 end
 
 
+facts("test solve elasticity increment, two elements") do
+
+    X = Float64[0 0; 1 0; 2 0; 0 1; 1 1; 2 1]'
+  elmap = [1 2 5 4; 2 3 6 5]'
+  nodalloads = [0 0; 0 0; 0 0; 0 0; 0 0; -3 0]'
+  @debug("nodal loads:\n", nodalloads)
+  dirichletbc = [0 0; NaN NaN; NaN NaN; 0 0; NaN NaN; NaN NaN]'
+
+    dim, nnodes = size(X)
+
+  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)
+
+    la = la*ones(1, nnodes)
+    mu = mu*ones(1, nnodes)
+    u = zeros(dim, nnodes)
+    du = zeros(dim, nnodes)
+
+  N(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
+    ]
+
+  dNdξ(ξ) = [-(1-ξ[2])/4.0    -(1-ξ[1])/4.0
+              (1-ξ[2])/4.0    -(1+ξ[1])/4.0
+              (1+ξ[2])/4.0     (1+ξ[1])/4.0
+             -(1+ξ[2])/4.0     (1-ξ[1])/4.0]
+
+  ipoints = 1/sqrt(3)*[-1 -1; 1 -1; 1 1; -1 1]
+  iweights = [1 1 1 1]
+
+  for i=1:10
+    solve_elasticity_increment!(X, u, du, elmap, nodalloads, dirichletbc,
+                                la, mu, N, dNdξ, ipoints, iweights)
+    @debug("increment:\n",du)
+    u += du
+    if norm(du) < 1.0e-9
+      break
+    end
+  end
+    @debug("solution\n",u)
+    # Known to fail, test against elmer.
+  @pending u[2, 6] => roughly(:something)
+end
+
+
 using JuliaFEM.elasticity_solver: interpolate
 facts("test interpolation of different field variables") do
     N(xi) = [