diff --git a/src/assembly.jl b/src/assembly.jl
index 703eaf5..39061b6 100644
--- a/src/assembly.jl
+++ b/src/assembly.jl
@@ -6,10 +6,10 @@
 type CAssembly
     interior_dofs :: Vector{Int}
     boundary_dofs :: Vector{Int}
-    F :: Factorization
+    F :: Union{Factorization, Matrix}
     Kc :: SparseMatrixCSC
     fc :: SparseMatrixCSC
-    Ki :: SparseMatrixCSC
+    Kib :: SparseMatrixCSC
     fi :: SparseMatrixCSC
 end
 
@@ -24,77 +24,98 @@ end
 
 function assemble(problem::AllProblems, time::Float64)
     assembly = Assembly()
-    for element in get_elements(problem)
+    ne = length(get_elements(problem))
+    p = ne > 10 ? round(Int, ne/10) : ne
+    for (i, element) in enumerate(get_elements(problem))
+        mod(i, p) == 0 && info("Assemble: ", round(Int, i/ne*100), " % done")
         assemble!(assembly, problem, element, time)
     end
     return assembly
 end
 
-""" Return condensed system. """
-function assemble(problem::FieldProblem, time::Float64, boundary_dofs::Vector{Int})
-    assembly = Assembly()
-    for element in get_elements(problem)
-        assemble!(assembly, problem, element, time)
-    end
-    return condensate(assembly, boundary_dofs)
-end
-
-function condensate(assembly::Assembly, boundary_dofs_::Vector{Int})
-    K = sparse(assembly.stiffness_matrix)
+""" Calculate reduced stiffness matrix. """
+function reduce(assembly::Assembly, boundary_dofs_::Vector{Int})
     all_dofs = unique(assembly.stiffness_matrix.I)
     boundary_dofs = intersect(all_dofs, boundary_dofs_)
     interior_dofs = setdiff(all_dofs, boundary_dofs_)
+
+    K = sparse(assembly.stiffness_matrix)
+    f = sparse(assembly.force_vector)
+
     dim = size(K, 1)
-    f = sparse(assembly.force_vector, dim, 1)
+
+    # empty assembly to release memory for factorization
+    empty!(assembly.stiffness_matrix)
+    empty!(assembly.force_vector)
+
+    if dim < 100000
+        # no need to do any reduction of matrix size at all
+        return CAssembly([], all_dofs, Matrix{Float64}(), K, f, spzeros(0, 0), spzeros(0,1))
+    end
 
     # check that matrix is symmetric
-    asdf = maximum(abs(1/2*(K + K') - K))
-    if asdf > 1.0e-6
-        info(full(K))
-        error("asdf $asdf > 1.0e-6")
-    end
-    
+    s = maximum(abs(1/2*(K + K') - K))
+    @assert s < 1.0e-6
     K = 1/2*(K + K')
-    F::Factorization = cholfact(K[interior_dofs, interior_dofs])
 
-#   info("condensation: all dofs: ", all_dofs)
-#   info("condensation: interior dofs: ", interior_dofs)
-#   info("condensation: boundary dofs: ", boundary_dofs)
-#   info("manually condensated")
-#   Kman = K[boundary_dofs, boundary_dofs] - K[boundary_dofs,interior_dofs] * inv(full(K[interior_dofs, interior_dofs])) * K[interior_dofs, boundary_dofs]
-#   info("\n$(full(Kman))")
-    #info("K = \n$(full(K))")
-    #Ki = K[interior_dofs, boundary_dofs]
-    Ki = K[interior_dofs, boundary_dofs]
+    Kib = K[interior_dofs, boundary_dofs]
+    Kbb = K[boundary_dofs, boundary_dofs]
     fi = f[interior_dofs]
-#   info("condensated using factorization")
-#   LL = K[boundary_dofs, boundary_dofs] - K[boundary_dofs, interior_dofs] * (K[interior_dofs, interior_dofs] \ K[interior_dofs, boundary_dofs])
-#   info(LL)
+    fb = f[boundary_dofs]
 
-    Ks = F \ Ki
-    Fs = F \ fi
+    F = cholfact(K[interior_dofs, interior_dofs])
+    K = spzeros(0, 0)
 
-    dim = size(K, 1)
+#=
+    if dim < 100000
+        # for small problems we don't need to care about memory usage
+        Kd = Kib' * (F \ Kib)
+    else
+        # for larger problems calculate schur complement in pieces
+        nb = length(boundary_dofs)
+        p = nb > 10 ? round(Int, nb/10) : nb
+        Kd = zeros(nb, nb)
+        for bi in 1:nb
+            mod(bi, p) == 0 && info("Reduction: ", round(Int, bi/nb*100), " % done")
+            C = full(F \ Kib[:, bi])
+            for bj in 1:nb
+                d = Kib[:, bj]
+                Kd[bj,bi] = dot(C[rowvals(d)], nonzeros(d))
+            end
+        end
+    end
     Kc = spzeros(dim, dim)
+    Kc[boundary_dofs, boundary_dofs] = Kbb - Kd
+=#
+
+    chunks = round(Int, dim/3000)
+    info("Reduction is done in $chunks chunks.")
+    nb = length(boundary_dofs)
+    kk = round(Int, collect(linspace(0, nb, chunks+1)))
+    sl = [kk[j]+1:kk[j+1] for j=1:length(kk)-1]
+    Kc = spzeros(dim, dim)
+    for (k,sli) in enumerate(sl)
+        b1 = boundary_dofs[sli]
+        Sc = F \ Kib[:,sli]
+        for slj in sl
+            b2 = boundary_dofs[slj]
+            Kc[b2,b1] = Kbb[slj,sli] - Kib[:,slj]'*Sc
+        end
+        info("Reduction: ", round(k/chunks*100, 0), " % done")
+    end
+
     fc = spzeros(dim, 1)
-    Kc[boundary_dofs, boundary_dofs] = K[boundary_dofs, boundary_dofs] - Ki' * Ks
-    fc[boundary_dofs] = f[boundary_dofs] - Ki' * Fs
+    fc[boundary_dofs] = fb - Kib' * (F \ fi)
 
-
-    return CAssembly(interior_dofs, boundary_dofs, F, Kc, fc, Ki, fi)
+    return CAssembly(interior_dofs, boundary_dofs, F, Kc, fc, Kib, fi)
 end
 
 function reconstruct!(ca::CAssembly, x::SparseMatrixCSC)
-#   info("size of la = ", size(la))
-#   info("size of ca.Ki = ", size(ca.Ki))
-#   info("size of ca.fi = ", size(ca.fi))
-#   info("size of la[ca.interior_dofs] = ", size(la[ca.interior_dofs]))
-#   info("interior dofs: $(ca.interior_dofs)")
-#   info("boundary dofs: $(ca.boundary_dofs)")
-#   info("ca.fi = $(ca.fi')")
-#   info("sol1 = ", full(ca.F \ ca.fi)')
-#   info("sol2 = ", full(ca.F \ (ca.Ki*x[ca.boundary_dofs]))')
-    x[ca.interior_dofs] += ca.F \ (ca.fi - ca.Ki*x[ca.boundary_dofs])
+    if isa(ca.F, Factorization)
+        x[ca.interior_dofs] = ca.F \ (ca.fi - ca.Kib*x[ca.boundary_dofs])
+    else  # normal inverse of matrix
+        x[ca.interior_dofs] = ca.F * (ca.fi - ca.Kib*x[ca.boundary_dofs])
+    end
 end
 
 function Base.(:+)(ass1::Assembly, ass2::Assembly)
diff --git a/src/core.jl b/src/core.jl
index 427a304..0cf93ce 100644
--- a/src/core.jl
+++ b/src/core.jl
@@ -87,3 +87,4 @@ include("directsolver.jl") # parallel sparse direct solver for non-linear proble
 ### MORTAR STUFF ###
 include("mortar.jl")  # mortar projection
 
+include("abaqus_reader_old.jl")
diff --git a/src/directsolver.jl b/src/directsolver.jl
index 10f0334..e32838f 100644
--- a/src/directsolver.jl
+++ b/src/directsolver.jl
@@ -38,7 +38,7 @@ function time_elapsed(timing, what::ASCIIString)
 end
 
 """ Call solver to solve a set of problems. """
-function call(solver::DirectSolver, time::Number=0.0)
+function call(solver::DirectSolver, ::Type{Val{:noreduce}}, time::Number=0.0)
     #@assert length(solver.field_problems) == 1
     info("# of field problems: $(length(solver.field_problems))")
     info("# of boundary problems: $(length(solver.boundary_problems))")
@@ -229,3 +229,204 @@ function call(solver::DirectSolver, time::Number=0.0)
 
 end
 
+
+""" Call solver to solve a set of problems. """
+function call(solver::DirectSolver, time::Number=0.0)
+    info("# of field problems: $(length(solver.field_problems))")
+    info("# of boundary problems: $(length(solver.boundary_problems))")
+    @assert solver.nonlinear_problem == true
+    
+    timing = Dict{ASCIIString, Float64}()
+    tic(timing, "solver")
+    tic(timing, "initialization")
+
+    # check that all problems are "same kind"
+    field_name = get_unknown_field_name(solver.field_problems[1])
+    field_dim = get_unknown_field_dimension(solver.field_problems[1])
+    for field_problem in solver.field_problems
+        get_unknown_field_name(field_problem) == field_name || error("several different fields not supported yet")
+        get_unknown_field_dimension(field_problem) == field_dim || error("several different field dimensions not supported yet")
+    end
+
+    # create initial fields for this increment
+    # i.e., copy last known values as initial guess
+    # for this increment
+
+    for field_problem in solver.field_problems
+        for element in get_elements(field_problem)
+            gdofs = get_gdofs(element, field_dim)
+            if haskey(element, field_name)
+                if !isapprox(last(element[field_name]).time, time)
+                    last_data = copy(last(element[field_name]).data)
+                    push!(element[field_name], time => last_data)
+                end
+            else
+                data = Vector{Float64}[zeros(field_dim) for i in 1:length(element)]
+                element[field_name] = (time => data)
+            end
+        end
+    end
+
+    for boundary_problem in solver.boundary_problems
+        for element in get_elements(boundary_problem)
+            gdofs = get_gdofs(element, field_dim)
+            data = Vector{Float64}[zeros(field_dim) for i in 1:length(element)]
+            if haskey(element, "reaction force")
+                if !isapprox(last(element["reaction force"]).time, time)
+                    push!(element["reaction force"], time => data)
+                end
+            else
+                element["reaction force"] = (time => data)
+            end
+        end
+    end
+
+    toc(timing, "initialization")
+
+    dim = 0
+
+    for iter=1:solver.max_iterations
+        info("Starting iteration $iter")
+        tic(timing, "non-linear iteration")
+
+        mapper = solver.parallel ? pmap : map
+
+        info("Assembling boundary problems...")
+        tic(timing, "boundary assembly")
+        boundary_assembly = sum(mapper((p)->assemble(p, time), solver.boundary_problems))
+        boundary_dofs = unique(boundary_assembly.stiffness_matrix.I)
+        info("# of interface dofs: $(length(boundary_dofs))")
+        toc(timing, "boundary assembly")
+
+        info("Assembling field problems...")
+        dim = 0
+        assemblies = []
+        for (i, problem) in enumerate(solver.field_problems)
+            info("Assembling body $i...")
+            tic(timing, "field assembly")
+            field_assembly = assemble(problem, time)
+            toc(timing, "field assembly")
+            field_dofs = unique(field_assembly.stiffness_matrix.I)
+            info("# of dofs in problem $i: $(length(field_dofs))")
+            info("Eliminating interior dofs for body $i...")
+            dim = maximum([dim, maximum(field_dofs)])
+            tic(timing, "condensate")
+            cfield_assembly = reduce(field_assembly, boundary_dofs)
+            toc(timing, "condensate")
+            push!(assemblies, cfield_assembly)
+        end
+
+        tic(timing, "create sparse matrices")
+        K = spzeros(dim, dim)
+        f = spzeros(dim, 1)
+
+        for (i, assembly) in enumerate(assemblies)
+            resize!(assembly.Kc, dim, dim)
+            resize!(assembly.fc, dim, 1)
+            K += assembly.Kc
+            f += assembly.fc
+        end
+
+        C = sparse(boundary_assembly.stiffness_matrix, dim, dim)
+        g = sparse(boundary_assembly.force_vector, dim, 1)
+        A = [K C'; C spzeros(dim, dim)]
+        b = [f; g]
+        toc(timing, "create sparse matrices")
+
+        info("Solving interface system")
+        tic(timing, "solution of system")
+        nz = sort(unique(rowvals(A)))  # take only non-zero rows
+        sol = zeros(b)
+        sol[nz] = A[nz,nz] \ full(b[nz])
+        toc(timing, "solution of system")
+
+#=
+        try
+        catch
+            dump(round(full(A[nz,nz]), 3))
+            dump(round(full(b[nz]'), 3))
+            for (i, assembly) in enumerate(assemblies)
+                info("assembly $i dump")
+                dump(round(full(assembly.Kc), 3))
+                dump(round(full(assembly.fc), 3)')
+            end
+            info("matrix K")
+            dump(round(full(K), 3))
+            info("interface matrix")
+            dump(round(full(C), 3))
+            info("final assembly to solve:")
+            dump(round(full(A), 3))
+            dump(round(full(b'), 3))
+            info("nonzero dofs: $nz")
+            info("nonzero dofs removed:")
+            dump(round(full(A[nz,nz]), 3))
+            dump(round(full(b[nz]'), 3))
+            detsys = det(A[nz,nz])
+            info("determinant of system: $detsys")
+            error("Solving system failed.")
+        end
+=#
+
+        info("Solved, calculating interior dofs...")
+        tic(timing, "back substitute")
+        for assembly in assemblies
+            length(assembly.interior_dofs) != 0 || continue
+            reconstruct!(assembly, sol)
+        end
+        toc(timing, "back substitute")
+
+        la = sol[dim+1:end]
+        la = vec(full(la))
+        sol = vec(full(sol))
+
+        info("Problem solved. solution norm: $(norm(sol[1:dim]))")
+
+        tic(timing, "update element data")
+        # update elements in field problems
+        for field_problem in solver.field_problems
+            for element in get_elements(field_problem)
+                gdofs = get_gdofs(element, field_dim)
+                local_sol = sol[gdofs]  # incremental data for element
+                local_sol = reshape(local_sol, field_dim, length(element))
+                local_sol = Vector{Float64}[local_sol[:,i] for i=1:length(element)]
+                last(element[field_name]).data += local_sol  # <-- added
+            end
+        end
+
+        # update elements in boundary problems
+        for boundary_problem in solver.boundary_problems
+            for element in get_elements(boundary_problem)
+                gdofs = get_gdofs(element, field_dim)
+                local_sol = la[gdofs]
+                local_sol = reshape(local_sol, field_dim, length(element))
+                local_sol = Vector{Float64}[local_sol[:,i] for i=1:length(element)]
+                last(element["reaction force"]).data = local_sol  # <-- replaced
+            end
+        end
+        toc(timing, "update element data")
+        toc(timing, "non-linear iteration")
+
+        if true
+            info("timing info for non-linear iteration:")
+            info("boundary assembly       : ", time_elapsed(timing, "boundary assembly"))
+            info("field assembly          : ", time_elapsed(timing, "field assembly"))
+            info("reduce stiffness matrix : ", time_elapsed(timing, "condensate"))
+            info("create sparse matrices  : ", time_elapsed(timing, "create sparse matrices"))
+            info("solution of system      : ", time_elapsed(timing, "solution of system"))
+            info("update element data     : ", time_elapsed(timing, "update element data"))
+            info("non-linear iteration    : ", time_elapsed(timing, "non-linear iteration"))
+        end
+
+        if norm(sol[1:dim]) < solver.tol
+            toc(timing, "solver")
+            info("solver finished in ", time_elapsed(timing, "solver"), " seconds.")
+            return (iter, true)
+        end
+
+    end
+
+    info("Warning: did not coverge in $(solver.max_iterations) iterations!")
+    return (solver.max_iterations, false)
+
+end
+
diff --git a/src/elasticity.jl b/src/elasticity.jl
index 83a2591..0afdc0d 100644
--- a/src/elasticity.jl
+++ b/src/elasticity.jl
@@ -80,7 +80,7 @@ function get_residual_vector{P<:ElasticityProblem}(problem::Problem{P}, element:
         J = det(element, ip, time)
         T = J^-1*F*S*F'
         #ip["cauchy stress"] = T
-        ip["gl strain"] = E
+        #ip["gl strain"] = E
 
         r += F*S*dbasis
     end
diff --git a/test/test_directsolver.jl b/test/test_directsolver.jl
index b87e2a0..a2d68e8 100644
--- a/test/test_directsolver.jl
+++ b/test/test_directsolver.jl
@@ -56,7 +56,7 @@ function test_solver_multiple_dirichlet_bc()
     push!(solver, problem3)
 
     # launch solver
-    norm = solver(0.0)
+    #norm = solver(0.0)
     norm = solver(1.0)
     disp = e1("displacement", [1.0, 1.0], 1.0)
     info("displacement at tip: $disp")
@@ -174,6 +174,6 @@ function test_solver_multiple_bodies_multiple_dirichlet_bc()
 
 end
 
-# test_solver_multiple_bodies_multiple_dirichlet_bc()
+test_solver_multiple_bodies_multiple_dirichlet_bc()
 
 end