diff --git a/src/elasticity_solver.jl b/src/elasticity_solver.jl
index 7d4f8bc..da7898f 100644
--- a/src/elasticity_solver.jl
+++ b/src/elasticity_solver.jl
@@ -10,6 +10,56 @@ VERSION < v"0.4-" && using Docile
 # directly if needed or using general interface combining data model and
 # solver.
 
+
+@doc """
+Interpolate field variable using basis functions f for point ip.
+This function tries to be as general as possible and allows interpolating
+lot of different fields.
+
+Parameters
+----------
+field :: Array{Number, dim}
+  Field variable
+basis :: Function
+  Basis functions
+ip :: Array{Number, 1}
+  Point to interpolate
+""" ->
+function interpolate(field::Array{Float64,1}, basis::Function, ip)
+    result = dot(field, basis(ip))
+    return result
+end
+function interpolate(field::Array{Float64,2}, basis::Function, ip)
+    m, n = size(field)
+    bip = basis(ip)
+    tmp = size(bip)
+    if length(tmp) == 1
+        ndim = 1
+        nnodes = tmp[1]
+    else
+        ndim, nnodes = size(bip)
+    end
+    if ndim == 1
+        if n == nnodes
+            result = field * bip
+        elseif m == nnodes
+            result = field' * bip
+        end
+    else
+      if n == nnodes
+        result = bip' * field
+    elseif m == nnodes
+      result = bip' * field'
+    end
+    end
+    if length(result) == 1
+        result = result[1]
+    end
+    return result
+end
+
+
+
 @doc """
 Calculate local tangent stiffness matrix and residual force vector R = T - F
 """ ->
@@ -25,9 +75,13 @@ function calc_local_matrices!(X, u, R, Kt, N, dNdξ, λ_, μ_, ipoints, iweights
     w = iweights[m]
     ξ = ipoints[m, :]
     # interpolate material parameters from element node fields
-    λ = (λ_*N(ξ))[1]
-    μ = (μ_*N(ξ))[1]
-        Jᵀ = X*dNdξ(ξ)
+    #λ = (λ_*N(ξ))[1]
+    #μ = (μ_*N(ξ))[1]
+    #    Jᵀ = X*dNdξ(ξ)
+    #@debug("Jt:\n",Jᵀ)
+    λ = interpolate(λ_, N, ξ)
+    μ = interpolate(μ_, N, ξ)
+    Jᵀ = interpolate(X, dNdξ, ξ)'
         detJ = det(Jᵀ)
         ∇N = inv(Jᵀ)*dNdξ(ξ)'
         ∇u = u*∇N'
@@ -283,8 +337,7 @@ function solve_elasticity_increment!(X, u, du, elmap, nodalloads,
     # solution
     free_dofs = find(isnan(dirichletbc))
     #du[free_dofs] = Kt[free_dofs, free_dofs] \ -reshape(R, 8)[free_dofs]
-    # TODO: cholesky decomposition
-    du[free_dofs] = full(A) \ -full(b)
+    du[free_dofs] = lufact(A) \ -full(b)
 end
 
 end
diff --git a/test/test_elasticity_solver.jl b/test/test_elasticity_solver.jl
index 00a4c80..a7c6c8c 100644
--- a/test/test_elasticity_solver.jl
+++ b/test/test_elasticity_solver.jl
@@ -8,7 +8,7 @@ using JuliaFEM.elasticity_solver: solve_elasticity_increment!
 
 facts("test solve elasticity increment") do
 
-  X = [0 0; 10 0; 10 1; 0 1]'
+  X = [0.0 0.0; 10.0 0.0; 10.0 1.0; 0.0 1.0]'
   elmap = [1; 2; 3; 4]
   nodalloads = [0 0; 0 0; 0 -2; 0 0]'
   @debug("nodal loads:\n", nodalloads)
@@ -54,6 +54,34 @@ facts("test solve elasticity increment") do
 end
 
 
+using JuliaFEM.elasticity_solver: interpolate
+facts("test interpolation of different field variables") do
+    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]
+    F1 = [36.0, 36.0, 36.0, 36.0]
+    F2 = [36.0 36.0 36.0 36.0]
+    F3 = F2'
+    F4 = [0.0 0.0; 10.0 0.0; 10.0 1.0; 0.0 1.0]'
+    F5 = F4'
+
+    @fact interpolate(F1, N, [0.0, 0.0]) => 36.0
+    @fact interpolate(F2, N, [0.0, 0.0]) => 36.0
+    @fact interpolate(F3, N, [0.0, 0.0]) => 36.0
+    @fact interpolate(F4, N, [0.0, 0.0]) => [5.0; 0.5]
+    @fact interpolate(F5, N, [0.0, 0.0]) => [5.0; 0.5]
+    @fact interpolate(F5, dNdξ, [0.0, 0.0]) => [5.0 0.0; 0.0 0.5]
+end
+
+
+
 using JuliaFEM.elasticity_solver: assemble!
 
 facts("test assembly of global matrix for 1 dim/node case") do