diff --git a/src/elements.jl b/src/elements.jl
index c7e2ec7..759fe00 100644
--- a/src/elements.jl
+++ b/src/elements.jl
@@ -188,10 +188,10 @@ end
 
 """ Return dual basis transformation matrix Ae. """
 function get_dualbasis(element::Element, time::Real)
-    if length(element.dualbasis) == 0
+    if length(element.A) == 0
         nnodes = size(element, 2)
-        D = zeros(nnodes, nnodes)
-        M = zeros(nnodes, nnodes)
+        De = zeros(nnodes, nnodes)
+        Me = zeros(nnodes, nnodes)
         for ip in get_integration_points(element, Val{3})
             w = ip.weight
             J = get_jacobian(element, ip, time)
@@ -207,8 +207,8 @@ function get_dualbasis(element::Element, time::Real)
             De += w*diagm(vec(N))
             Me += w*N'*N
         end
-        element.D = D
-        element.M = M
+        element.D = De
+        element.M = Me
         element.A = De*inv(Me)
     end
     return element.D, element.M, element.A
diff --git a/src/problems.jl b/src/problems.jl
index 85ab797..6274185 100644
--- a/src/problems.jl
+++ b/src/problems.jl
@@ -271,6 +271,7 @@ end
 """ Find nodes corresponding to dofs. """
 function find_nodes_by_dofs(problem::Problem, dofs)
     dim = get_unknown_field_dimension(problem)
+    return find_nodes_by_dofs(dim, dofs)
 end
 function find_nodes_by_dofs(dim, dofs)
     nodes = Int64[]
diff --git a/src/solver_utils.jl b/src/solver_utils.jl
index bb97462..b98ff26 100644
--- a/src/solver_utils.jl
+++ b/src/solver_utils.jl
@@ -1,7 +1,8 @@
 # This file is a part of JuliaFEM.
 # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
 
-subscript(i) = map(repr(i)) do c
+function subscript(i)
+    map(repr(i)) do c
         c   ==  '1' ? '\u2081' :
         c   ==  '2' ? '\u2082' :
         c   ==  '3' ? '\u2083' :
@@ -12,8 +13,9 @@ subscript(i) = map(repr(i)) do c
         c   ==  '8' ? '\u2088' :
         c   ==  '9' ? '\u2089' :
         c   ==  '0' ? '\u2080' :
-        error("Unexpected Chatacter")
+        error("Unexpected character")
     end
+end
 
 function pretty_print_constraint_equation(a, b, g; char1="u", char2="λ")
     s = ""
@@ -95,17 +97,17 @@ function handle_overconstraint_error!(problem, nodes, all_dofs, C1_, C1, C2_, C2
     """
     function action1(node_id, dofs)
         dofs_ = get_related_dofs(intersect(dofs, all_dofs))
+        # this will fail with dofs > 2 for some yet unknown reason
+        length(dofs_) > 2 && return dofs_, false
+
         any(has_lagrange_coefficients(dofs_)) && return dofs_, false
         C = full([C2[dofs_, :]; C2_[dofs_, :]])
         d = full([g[dofs_]; g_[dofs_]])
-        # this will fail with dofs > 2 for some yet unknown reason
-        length(dofs_) > 2 && return dofs_, false
         info("rank = $(rank(C)), dofs = $(length(dofs_))")
         rank(C) != length(dofs_) && return dofs_, false
         C2_[dofs_,:] = C2[dofs_,:] = 0
         g_[dofs_,:] = g[dofs_,:] = 0
         x = C \ d
-        x[abs(x) .< 1.0e-12] = 0
         C2[dofs_, dofs_] = eye(length(dofs_))
         g[dofs_] = x
         return dofs_, true
diff --git a/src/solvers.jl b/src/solvers.jl
index 80fcb9b..fe65c90 100644
--- a/src/solvers.jl
+++ b/src/solvers.jl
@@ -138,11 +138,12 @@ end
 # Tuple{Symbol,Any,Any} or Function
 
 type Solver
-    name :: ASCIIString
-    time :: Real
-    iteration :: Int
+    name :: ASCIIString # some descriptive name for problem
+    time :: Real        # current time
+    iteration :: Int    # iteration counter
+    ndofs :: Int        # total dimension of global stiffness matrix, i.e., dim*nnodes
     problems :: Vector{Problem}
-    is_linear_system :: Bool
+    is_linear_system :: Bool # setting this to true makes assumption of one step convergence
     nonlinear_system_max_iterations :: Int64
     nonlinear_system_convergence_tolerance :: Float64
     linear_system_solver :: Symbol
@@ -150,11 +151,12 @@ end
 
 function Solver(name::ASCIIString="default solver", time::Real=0.0)
     return Solver(
-        name,  # name
-        time,  # time
-        0,     # iteration counter
+        name,
+        time,
+        0,     # iteration #
+        0,     # ndofs
         [],    # array of problems
-        false, # is this a linear system which can be solved in a single iteration?
+        false, # is_linear_system
         10,    # max nonlinear iterations
         5.0e-5, # nonlinear iteration convergence tolerance
         :DirectLinearSolver # linear system solution method
@@ -210,6 +212,7 @@ function get_mortar_problems(solver::Solver)
     filter(is_mortar_problem, solver.problems)
 end
 
+
 """Return one combined field assembly for a set of field problems.
 
 Parameters
@@ -218,7 +221,7 @@ solver :: Solver
 
 Returns
 -------
-K, f :: SparseMatrixCOO
+K, f :: SparseMatrixCSC
 
 Notes
 -----
@@ -227,42 +230,65 @@ problems must have unique node ids.
 
 """
 function get_field_assembly(solver::Solver)
-    return get_field_assembly(get_field_problems(solver))
-end
-function get_field_assembly(problems::Vector{Problem})
+    problems = get_field_problems(solver)
     K = SparseMatrixCOO()
     f = SparseMatrixCOO()
     for problem in problems
         append!(K, problem.assembly.K)
         append!(f, problem.assembly.f)
     end
+    K = sparse(K)
+    solver.ndofs = size(K, 1)
+    f = sparse(f, solver.ndofs, 1)
     return K, f
 end
 
+
 """ Return one combined boundary assembly for a set of boundary problems.
 
 Returns
 -------
-C1, C2, D, g :: SparseMatrixCOO
+C1, C2, D, g :: SparseMatrixCSC
+
+Notes
+-----
+When some dof is constrained by multiple boundary problems an algorithm is
+launched what tries to do it's best to solve issue. It's far from perfect
+but is able to handle some basic situations occurring in corner nodes and
+crosspoints.
 
 """
 function get_boundary_assembly(solver::Solver)
-    return get_boundary_assembly(get_boundary_problems(solver))
-end
-function get_boundary_assembly(problems::Vector{Problem})
-    C1 = SparseMatrixCOO()
-    C2 = SparseMatrixCOO()
-    D = SparseMatrixCOO()
-    g = SparseMatrixCOO()
-    for problem in problems
-        append!(C1, problem.assembly.C1)
-        append!(C2, problem.assembly.C2)
-        append!(D, problem.assembly.D)
-        append!(g, problem.assembly.g)
+    ndofs = solver.ndofs
+    @assert ndofs != 0
+    C1 = spzeros(ndofs, ndofs)
+    C2 = spzeros(ndofs, ndofs)
+    D = spzeros(ndofs, ndofs)
+    g = spzeros(ndofs, 1)
+    for problem in get_boundary_problems(solver)
+        assembly = problem.assembly
+        C1_ = sparse(assembly.C1, ndofs, ndofs)
+        C2_ = sparse(assembly.C2, ndofs, ndofs)
+        D_ = sparse(assembly.D, ndofs, ndofs)
+        g_ = sparse(assembly.g, ndofs, 1)
+        already_constrained = get_nonzero_rows(C2)
+        new_constraints = get_nonzero_rows(C2_)
+        overconstrained_dofs = intersect(already_constrained, new_constraints)
+        if length(overconstrained_dofs) != 0
+            overconstrained_dofs = sort(overconstrained_dofs)
+            overconstrained_nodes = find_nodes_by_dofs(problem, overconstrained_dofs)
+            handle_overconstraint_error!(problem, overconstrained_nodes,
+                overconstrained_dofs, C1, C1_, C2, C2_, D, D_, g, g_)
+        end
+        C1 += C1_
+        C2 += C2_
+        D += D_
+        g += g_
     end
     return C1, C2, D, g
 end
 
+
 """ Solve linear system using LU factorization (UMFPACK).
 """
 function solve_linear_system(solver::Solver, ::Type{Val{:DirectLinearSolver}})
@@ -271,31 +297,26 @@ function solve_linear_system(solver::Solver, ::Type{Val{:DirectLinearSolver}})
 
     # assemble field problems
     K, f = get_field_assembly(solver)
-    K = sparse(K)
-    dim = size(K, 1)
-    f = sparse(f, dim, 1)
 
     # assemble boundary problems
     C1, C2, D, g = get_boundary_assembly(solver)
-    C1 = sparse(C1, dim, dim)
-    C2 = sparse(C2, dim, dim)
-    D = sparse(D, dim, dim)
-    g = sparse(g, dim, 1)
 
     # construct global system Ax=b and solve using lu factorization
     A = [K C1'; C2 D]
     b = [f; g]
-    nz1 = sort(unique(rowvals(A)))
-    nz2 = sort(unique(rowvals(A')))
-    x = zeros(length(b))
-    x[nz1] = lufact(A[nz1,nz2]) \ full(b[nz1])
 
-    u = x[1:dim]
-    la = x[dim+1:end]
+    nz = get_nonzero_rows(A)
+    x = zeros(length(b))
+    x[nz] = lufact(A[nz,nz]) \ full(b[nz])
+
+    ndofs = solver.ndofs
+    u = x[1:ndofs]
+    la = x[ndofs+1:end]
     info("UMFPACK: solved in ", time()-t0, " seconds. norm = ", norm(u))
     return u, la
 end
 
+
 """ Check convergence of problems.
 
 Notes
diff --git a/src/sparse.jl b/src/sparse.jl
index 0612739..4445258 100644
--- a/src/sparse.jl
+++ b/src/sparse.jl
@@ -135,7 +135,7 @@ Returns
 
 Ordered list of row indices.
 """
-function get_nonzero_rows(A::SparseMatrixCOO)
+function get_nonzero_rows(A::SparseMatrixCSC)
     # FIXME: This is probably a very inefficient way to do this.
     return sort(unique(rowvals(A)))
 end
diff --git a/test/test_node_dof_mapping.jl b/test/test_node_dof_mapping.jl
index 4b5c06b..1a9c0d2 100644
--- a/test/test_node_dof_mapping.jl
+++ b/test/test_node_dof_mapping.jl
@@ -17,5 +17,11 @@ end
     dim = 3
     nodes = find_nodes_by_dofs(dim, dofs)
     @test nodes == [1, 3]
+
+    dofs = [2, 12]
+    dim = 2
+    nodes = find_nodes_by_dofs(dim,dofs)
+    @test nodes == [1, 6]
+
 end