diff --git a/src/elasticity_solver.jl b/src/elasticity_solver.jl
index 7173636..3026bba 100644
--- a/src/elasticity_solver.jl
+++ b/src/elasticity_solver.jl
@@ -135,6 +135,89 @@ function assemble!(fe, eldofs_, I, V)
 end
 
 
+@doc """
+Eliminate Dirichlet boundary conditions from matrix
+
+Parameters
+----------
+dirichletbc : array [dim x nnodes]
+I, J, V : sparse matrix arrays
+
+Returns
+-------
+I, J, V : boundary conditions removed
+
+Notes
+-----
+pros:
+- matrix assembly remains positive definite
+cons:
+- maybe inefficient because of extra sparse matrix operations. (It's hard to remove stuff from sparse matrix.)
+- if u != 0 in dirichlet boundary requires extra care
+
+Raises
+------
+Exception, if displacement boundary conditions given, i.e.
+DX=2 for some node, for example.
+
+"""
+function eliminate_boundary_conditions(dirichletbc, I, J, V)
+    if any(dirichletbc .> 0)
+        throw("displacement boundary condition not supported")
+    end
+    # dofs to remove
+    free_dofs = find(isnan(dirichletbc))
+    remove_dofs = find(!isnan(dirichletbc))
+    @debug("Removing dofs: ", remove_dofs)
+    # this can be done more clever by removing corresponging indexes from I, J, and V
+    A = sparse(I, J, V)
+    A = A[free_dofs, free_dofs]
+    return findnz(A)
+end
+
+@doc """
+Eliminate Dirichlet boundary conditions from vector
+
+Parameters
+----------
+dirichletbc : array [dim x nnodes]
+I, V : sparse vector arrays
+
+Returns
+-------
+I, V : boundary conditions removed
+
+Notes
+-----
+pros:
+- matrix assembly remains positive definite
+cons:
+- maybe inefficient because of extra sparse matrix operations. (It's hard to remove stuff from sparse matrix.)
+- if u != 0 in dirichlet boundary requires extra care
+
+Raises
+------
+Exception, if displacement boundary conditions given, i.e.
+DX=2 for some node, for example.
+""" ->
+function eliminate_boundary_conditions(dirichletbc, I, V)
+    if any(dirichletbc .> 0)
+        throw("displacement boundary condition not supported")
+    end
+    # dofs to remove
+    free_dofs = find(isnan(dirichletbc))
+    remove_dofs = find(!isnan(dirichletbc))
+    @debug("Removing dofs: ", remove_dofs)
+    # this can be done more clever by removing corresponging indexes from I, J, and V
+    A = sparsevec(I, V)
+    A = A[free_dofs]
+    @debug("new vector: ", A)
+    i, j, v = findnz(A)
+    return i, v
+end
+
+
+
 @doc """
 Solve one increment of elasticity problem
 """ ->
diff --git a/test/test_elasticity_solver.jl b/test/test_elasticity_solver.jl
index 9c508ec..bffe381 100644
--- a/test/test_elasticity_solver.jl
+++ b/test/test_elasticity_solver.jl
@@ -119,3 +119,38 @@ facts("test assembly of global vector for 2 dim/node case") do
     expected = [1.0 2.0 5.0 8.0 6.0 8.0]'
     @fact S => expected
 end
+
+
+using JuliaFEM.elasticity_solver: eliminate_boundary_conditions
+facts("remove boundary conditions from matrix with 2 dof/node") do
+    # create sparse matrix 4x4 with some data
+    # 4x4 Array{Int64,2}:
+    #  1  5   9  13
+    #  2  6  10  14
+    #  3  7  11  15
+    #  4  8  12  16
+    A = sparse(reshape(1:4*4, 4, 4))
+    I, J, V = findnz(A)
+    # we plan to eliminate first dof of first node and second dof of second node
+    # expected output would be
+    # 6 10
+    # 7 11
+    dirichletbc = [0 NaN; NaN 0]'
+    I, J, V = eliminate_boundary_conditions(dirichletbc, I, J, V)
+    A2 = full(sparse(I, J, V))
+    @fact A2 => [6 10; 7 11]
+end
+
+facts("remove boundary conditions from vector with 2 dof/node") do
+    # create sparse vector dim 4 with some data
+    #  1 2 3 4 '
+    A = sparsevec([1, 2, 3, 4])
+    I, J, V = findnz(A)
+    # we plan to eliminate first dof of first node and second dof of second node
+    # expected output would be
+    # 2 3
+    dirichletbc = [0 NaN; NaN 0]'
+    I, V = eliminate_boundary_conditions(dirichletbc, I, V)
+    A2 = full(sparsevec(I, V))
+    @fact A2 => [2 3]'
+end