diff --git a/REQUIRE b/REQUIRE
index 251912c..288935f 100644
--- a/REQUIRE
+++ b/REQUIRE
@@ -7,3 +7,4 @@ JLD
 Compat
 DataFrames
 Formatting
+Logging
diff --git a/src/JuliaFEM.jl b/src/JuliaFEM.jl
index 416f2a3..008f74e 100644
--- a/src/JuliaFEM.jl
+++ b/src/JuliaFEM.jl
@@ -1,6 +1,12 @@
 # This file is a part of JuliaFEM.
 # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
 
+using Logging
+
+if haskey(ENV, "JULIAFEM_LOGLEVEL")
+    ENV["JULIAFEM_LOGLEVEL"] == "DEBUG" && Logging.configure(level=DEBUG)
+end
+
 """
 This is JuliaFEM -- Finite Element Package
 """
diff --git a/src/elements.jl b/src/elements.jl
index 17c46dd..98d43ab 100644
--- a/src/elements.jl
+++ b/src/elements.jl
@@ -11,10 +11,12 @@ type Element{E<:AbstractElement}
     properties :: E
 end
 
-function Element{E<:AbstractElement}(::Type{E}, connectivity=[], integration_points=[], id=-1, fields=Dict(), properties...)
-    variant = E(properties...)
-    element = Element{E}(id, connectivity, integration_points, fields, variant)
-    return element
+function Element{E<:AbstractElement}(::Type{E}, id::Int64, connectivity::Vector{Int64})
+    return Element{E}(id, connectivity, [], Dict(), E())
+end
+
+function Element{E<:AbstractElement}(::Type{E}, connectivity::Vector{Int64})
+    return Element{E}(-1, connectivity, [], Dict(), E())
 end
 
 function getindex(element::Element, field_name::AbstractString)
@@ -59,9 +61,15 @@ function call(element::Element, ip, time::Float64=0.0)
 end
 
 function call(element::Element, ip, time::Float64, ::Type{Val{:Jacobian}})
-    X = element["geometry"](time)
+    X = element("geometry", time)
     dN = get_dbasis(element, ip, time)
-    J = sum([kron(dN[:,i], X[i]') for i=1:length(X)])
+    nbasis = length(element)
+    if isa(X.data, Vector)
+        J = sum([kron(dN[:,i], X[i]') for i=1:nbasis])
+    else
+        c = get_connectivity(element)
+        J = sum([kron(dN[:,i], X[c[i]]') for i=1:nbasis])
+    end
     return J
 end
 
@@ -122,11 +130,16 @@ function call(element::Element, field::Field, ip, time::Float64)
     field_ = field(time)
     basis = element(ip, time)
     n = length(element)
-    m = length(field_)
-    if n != m
-        error("Error when trying to interpolate field $field at coords $ip and time $time: element length is $n and field length is $m, f = Nᵢfᵢ makes no sense!")
+    if isa(field_.data, Vector)
+        m = length(field_)
+        if n != m
+            error("Error when trying to interpolate field $field at coords $ip and time $time: element length is $n and field length is $m, f = Nᵢfᵢ makes no sense!")
+        end
+        return sum([field_[i]*basis[i] for i=1:n])
+    else
+        c = get_connectivity(element)
+        return sum([field_[c[i]]*basis[i] for i=1:n])
     end
-    return sum([field_[i]*basis[i] for i=1:n])
 end
 
 function size(element::Element, dim)
@@ -145,13 +158,19 @@ As a result element now have time invariant (variable) vector field "geometry" w
 
 """
 function update!(element::Element, field_name, data::Dict)
-    element[field_name] = [data[i] for i in get_connectivity(element)]
+    element[field_name] = Field(data)
+    #element[field_name] = [data[i] for i in get_connectivity(element)]
 end
 
 function update!{K,V}(element::Element, field_name, data::Pair{Float64, Dict{K, V}})
-    time, field_data = data
-    element_data = V[field_data[i] for i in get_connectivity(element)]
-    update!(element, field_name, time => element_data)
+    #time, field_data = data
+    #element_data = V[field_data[i] for i in get_connectivity(element)]
+    #update!(element, field_name, time => element_data)
+    if haskey(element, field_name)
+        update!(element[field_name], data)
+    else
+        element[field_name] = Field(data)
+    end
 end
 
 function update!(element::Element, field_name::AbstractString, datas::Union{Real, Vector, Pair{Float64, Union{Float64, Real, Vector{Any}}}}...)
@@ -321,4 +340,3 @@ function inside{E}(element::Element{E}, X, time)
     xi = get_local_coordinates(element, X, time)
     return inside(E, xi)
 end
-
diff --git a/src/fields.jl b/src/fields.jl
index f8ed44d..4425128 100644
--- a/src/fields.jl
+++ b/src/fields.jl
@@ -97,11 +97,24 @@ end
 function Field{T}(data::Pair{Float64, T}...)
     return DCTV([Increment{T}(d[1], d[2]) for d in data])
 end
-
+#=
 function Field{T}(data::Pair{Float64, Vector{T}}...)
     return DVTV([Increment{Vector{T}}(d[1], d[2]) for d in data])
 end
 
+function Field{T}(data::Pair{Float64, Dict{Int64, T}}...)
+    return DVTV([Increment{Dict{Int64, T}}(d[1], d[2]) for d in data])
+end
+=#
+
+function Field{T<:Union{Vector, Dict}}(data::Pair{Float64, T}...)
+    return DVTV([Increment{T}(d[1], d[2]) for d in data])
+end
+
+function Field(data::Dict)
+    return DVTI(data)
+end
+
 function convert{T}(::Type{DCTV}, data::Pair{Real, Vector{T}}...)
     return DCTV([Increment{Vector{T}}(d[1], d[2]) for d in data])
 end
@@ -222,11 +235,11 @@ function Base.(:*)(N::Matrix, f::DCTI)
 end
 
 
-#  
+#
 #   Multiply DVTI field with another vector T. Vector length
 #   must match to the field length and this can be used mainly
 #   for interpolation purposes, i.e., u = ∑ Nᵢuᵢ
-#  
+#
 function Base.(:*)(T::Vector, f::DVTI)
     @assert length(T) == length(f)
     return sum([T[i]*f[i] for i=1:length(f)])
@@ -430,3 +443,7 @@ end
 function keys(field::DVTV)
     return Float64[increment.time for increment in field]
 end
+
+function setindex!(field::Field, val, idx::Int64)
+    field.data[idx] = val
+end
diff --git a/src/io.jl b/src/io.jl
index 9a07c0c..e90c808 100644
--- a/src/io.jl
+++ b/src/io.jl
@@ -22,6 +22,42 @@ function haskey(x::XMLElement, key::AbstractString)
     return has_child(x, key) || has_attribute(x, key)
 end
 
+function has_child(x::XMLElement, child_name::AbstractString)
+    return get_child(x, child_name) != nothing
+end
+
+function get_attribute(x::XMLElement, attr_name::AbstractString)
+    attr = attribute(x, attr_name)
+    numeric = tryparse(Int64, attr)
+    isnull(numeric) && (numeric = tryparse(Float64, attr))
+    isnull(numeric) && return attr
+    return get(numeric)
+end
+
+function new_child(xparent::XMLElement, name::AbstractString, attrs::Dict)
+    x = new_child(xparent, name)
+    for (k, v) in attrs
+        x[k] = v
+    end
+    return x
+end
+
+function new_child(xparent::XMLElement, name::AbstractString, attrs::Pair...)
+    x = new_child(xparent, name)
+    for (k, v) in attrs
+        x[k] = v
+    end
+    return x
+end
+
+""" Basic traverse support, so that it's possible to find data from xml using
+path syntax e.g. /foo/bar[2]/baz[@Name=Frame 1]/DataItem. If several elements
+with same name exists in tree, pick first by default and next ones can be picked
+using [] syntax or [@attr=value] syntax, see [1] for details. For last item use
+[end].
+
+[1] http://www.xdmf.org/index.php/XDMF_Model_and_Format
+"""
 function get_child(x::XMLElement, child_name::AbstractString)
     '/' in child_name && return nothing
     m = match(r"(\w+)\[(.+)\]", child_name)
@@ -52,18 +88,6 @@ function get_child(x::XMLElement, child_name::AbstractString)
     throw("Unable to parse: $(m[2])")
 end
 
-function has_child(x::XMLElement, child_name::AbstractString)
-    return get_child(x, child_name) != nothing
-end
-
-function get_attribute(x::XMLElement, attr_name::AbstractString)
-    attr = attribute(x, attr_name)
-    numeric = tryparse(Int64, attr)
-    isnull(numeric) && (numeric = tryparse(Float64, attr))
-    isnull(numeric) && return attr
-    return get(numeric)
-end
-
 function getindex(x::XMLElement, attr_name::AbstractString)
     attr_name = strip(attr_name, '/')
     child = get_child(x, attr_name)
@@ -82,22 +106,6 @@ function getindex(x::XMLElement, attr_name::AbstractString)
     end
 end
 
-function new_child(xparent::XMLElement, name::AbstractString, attrs::Dict)
-    x = new_child(xparent, name)
-    for (k, v) in attrs
-        x[k] = v
-    end
-    return x
-end
-
-function new_child(xparent::XMLElement, name::AbstractString, attrs::Pair...)
-    x = new_child(xparent, name)
-    for (k, v) in attrs
-        x[k] = v
-    end
-    return x
-end
-
 type Xdmf
     name :: AbstractString
     xml :: XMLElement
diff --git a/src/postprocess_utils.jl b/src/postprocess_utils.jl
index a861c62..b6a24e7 100644
--- a/src/postprocess_utils.jl
+++ b/src/postprocess_utils.jl
@@ -199,28 +199,6 @@ function copy_field!(src_problem::Problem, dst_problem::Problem, field_name, tim
     copy_field!(src_problem.elements, dst_problem.elements, field_name, time)
 end
 
-""" Return field calculated to nodal points for elements in problem p. """
-function call(problem::Problem, field_name::AbstractString, time::Float64=0.0)
-    f = nothing
-    for element in get_elements(problem)
-        haskey(element, field_name) || continue
-        for (c, v) in zip(get_connectivity(element), element(field_name, time))
-            if f == nothing
-                f = Dict(c => v)
-                continue
-            end
-            if haskey(f, c)
-                if !isapprox(f[c], v)
-                    info("several values for single node when returning field $field_name")
-                    info("already have: $(f[c]), and trying to set $v")
-                end
-            else
-                f[c] = v
-            end
-        end
-    end
-    return f
-end
 
 function to_dataframe(u::Dict, abbreviation::Symbol)
     length(u) != 0 || return DataFrame()
@@ -419,20 +397,3 @@ function calculate_second_moment_of_mass(problem::Problem, X=[0.0, 0.0, 0.0], ti
     end
     return I
 end
-
-function getindex(problem::Problem, field_name::AbstractString)
-    info("fetching result $field_name")
-    timeframes = []
-    for frame in first(problem.elements)[field_name].data
-        push!(timeframes, frame.time)
-    end
-    info("time frames: $timeframes")
-    conn = get_connectivity(problem)
-    increments = Increment[]
-    for time in timeframes
-        p = problem(field_name, time)
-        data = [p[id] for id in conn]
-        push!(increments, Increment(time, data))
-    end
-    return DVTV(increments)
-end
diff --git a/src/problems.jl b/src/problems.jl
index 9514651..b3b0ff9 100644
--- a/src/problems.jl
+++ b/src/problems.jl
@@ -8,8 +8,8 @@ abstract MixedProblem <: AbstractProblem
 
 """
 General linearized problem to solve
-    (K₁+K₂)*Δu +   C1.T*λ = f₁+f₂
-         C2*Δu +      D*λ = g
+    (K₁+K₂)Δu  +   C1*Δλ = f₁+f₂
+         C2Δu  +    D*Δλ = g
 """
 type Assembly
 
@@ -19,7 +19,7 @@ type Assembly
     K :: SparseMatrixCOO   # stiffness matrix
     Kg :: SparseMatrixCOO  # geometric stiffness matrix
     f :: SparseMatrixCOO   # force vector
-    fg :: SparseMatrixCOO  # 
+    fg :: SparseMatrixCOO  #
 
     # for boundary assembly
     C1 :: SparseMatrixCOO
@@ -90,6 +90,7 @@ type Problem{P<:AbstractProblem}
     elements :: Vector{Element}
     dofmap :: Dict{Element, Vector{Int64}} # connects element local dofs to global dofs
     assembly :: Assembly
+    fields :: Dict{AbstractString, Field}
     properties :: P
 end
 
@@ -104,10 +105,10 @@ julia> prob2 = Problem(Elasticity, 3)
 
 """
 function Problem{P<:FieldProblem}(::Type{P}, name::AbstractString, dimension::Int64)
-    return Problem{P}(name, dimension, "none", [], Dict(), Assembly(), P())
+    return Problem{P}(name, dimension, "none", [], Dict(), Assembly(), Dict(), P())
 end
 function Problem{P<:FieldProblem}(::Type{P}, dimension::Int64)
-    return Problem{P}("$P problem", dimension, "none", [], Dict(), Assembly(), P())
+    return Problem{P}("$P problem", dimension, "none", [], Dict(), Assembly(), Dict(), P())
 end
 
 """ Construct a new boundary problem.
@@ -117,16 +118,16 @@ Examples
 Create Dirichlet boundary problem for vector-valued (dim=3) elasticity problem.
 
 julia> bc1 = Problem(Dirichlet, "support", 3, "displacement")
-
+solver.
 """
 function Problem{P<:BoundaryProblem}(::Type{P}, name, dimension, parent_field_name)
-    return Problem{P}(name, dimension, parent_field_name, [], Dict(), Assembly(), P())
+    return Problem{P}(name, dimension, parent_field_name, [], Dict(), Assembly(), Dict(), P())
 end
 function Problem{P<:BoundaryProblem}(::Type{P}, main_problem::Problem)
     name = "$P problem"
     dimension = get_unknown_field_dimension(main_problem)
     parent_field_name = get_unknown_field_name(main_problem)
-    return Problem{P}(name, dimension, parent_field_name, [], Dict(), Assembly(), P())
+    return Problem{P}(name, dimension, parent_field_name, [], Dict(), Assembly(), Dict(), P())
 end
 
 function get_formulation_type{P<:FieldProblem}(problem::Problem{P})
@@ -198,6 +199,7 @@ function update!(problem::Problem, assembly::Assembly, u::Vector, la::Vector; ve
     # incremental formulation we solve KΔu = f and u = u + Δu
     assembly.u_prev = copy(assembly.u)
     assembly.la_prev = copy(assembly.la)
+
     if get_formulation_type(problem) == :total
         verbose && info("$(problem.name): total formulation, replacing solution vector with new values")
         assembly.u = u
@@ -225,7 +227,7 @@ end
 
 Notes
 -----
-If length of solution vector != number of nodes, i.e. field dimension is 
+If length of solution vector != number of nodes, i.e. field dimension is
 something other than 1, reshape vectors so it's length matches to the
 number of nodes so that one can easily get nodal results.
 """
@@ -282,9 +284,50 @@ function length(problem::Problem)
 end
 
 function update!(problem::Problem, field_name::AbstractString, data)
+    if haskey(problem.fields, field_name)
+        update!(problem.fields[field_name], field_name::AbstractString, data)
+    else
+        problem.fields[field_name] = Field(data)
+    end
     update!(problem.elements, field_name::AbstractString, data)
 end
 
+function haskey(problem::Problem, field_name::AbstractString)
+    return haskey(problem.fields, field_name)
+end
+
+function getindex(problem::Problem, field_name::AbstractString)
+    return problem.fields[field_name]
+end
+
+""" Return field calculated to nodal points for elements in problem p. """
+function call(problem::Problem, field_name::AbstractString, time::Float64=0.0)
+    if haskey(problem, field_name)
+        return problem[field_name](time)
+    end
+    f = nothing
+    for element in get_elements(problem)
+        haskey(element, field_name) || continue
+        for (c, v) in zip(get_connectivity(element), element(field_name, time))
+            if f == nothing
+                f = Dict(c => v)
+                continue
+            end
+            if haskey(f, c)
+                if !isapprox(f[c], v)
+                    info("several values for single node when returning field $field_name")
+                    info("already have: $(f[c]), and trying to set $v")
+                end
+            else
+                f[c] = v
+            end
+        end
+    end
+    f == nothing && return f
+    update!(problem, field_name, time => f)
+    return f
+end
+
 """ Return the dimension of the unknown field of this problem. """
 function get_unknown_field_dimension(problem::Problem)
     return problem.dimension
@@ -381,4 +424,3 @@ function find_nodes_by_dofs(dim, dofs)
     end
     return nodes
 end
-
diff --git a/src/problems_elasticity.jl b/src/problems_elasticity.jl
index 81c5a71..b7dc871 100644
--- a/src/problems_elasticity.jl
+++ b/src/problems_elasticity.jl
@@ -4,7 +4,7 @@
 """ Elasticity equations.
 
 Field equation is:
-    
+
     m∂²u/∂t² = ∇⋅σ - b
 
 Weak form is: find u∈U such that ∀v in V
@@ -151,7 +151,7 @@ function assemble{El<:Elasticity2DVolumeElements}(problem::Problem{Elasticity},
         # cauchy_stress = F'*stress*F/det(F)
         # cauchy_stress = [cauchy_stress[1,1]; cauchy_stress[2,2]; cauchy_stress[1,2]]
         # update!(ip, "cauchy stress", time => cauchy_stress)
-        
+
         # material stiffness end
 
         if props.geometric_stiffness
@@ -170,13 +170,13 @@ function assemble{El<:Elasticity2DVolumeElements}(problem::Problem{Elasticity},
             S2[2,2] = stress_vec[2]
             S2[1,2] = S2[2,1] = stress_vec[3]
             S2[3:4,3:4] = S2[1:2,1:2]
-            
+
             Kg += w*BNL'*S2*BNL # geometric stiffness
 
         end
 
         # rhs, internal and external load
-        
+
         f -= w*BL'*stress_vec
 
         if haskey(element, "displacement load")
@@ -747,4 +747,3 @@ function call(problem::Problem, element::Element, ip, time::Float64, ::Type{Val{
     haskey(element, "geometry") || return nothing
     return element("geometry", ip, time)
 end
-
diff --git a/src/problems_heat.jl b/src/problems_heat.jl
index 2716b55..c7956b0 100644
--- a/src/problems_heat.jl
+++ b/src/problems_heat.jl
@@ -253,4 +253,3 @@ function assemble!{E<:Heat2DSurfaceElements}(assembly::Assembly, problem::Proble
     add!(assembly.K, gdofs, gdofs, K)
     add!(assembly.f, gdofs, fq)
 end
-
diff --git a/src/problems_mortar_2d.jl b/src/problems_mortar_2d.jl
index de42cdf..7934bd4 100644
--- a/src/problems_mortar_2d.jl
+++ b/src/problems_mortar_2d.jl
@@ -175,7 +175,7 @@ function assemble!(problem::Problem{Mortar}, time::Float64, ::Type{Val{1}}, ::Ty
                 n_s = N1*n1 # normal direction in gauss point
                 xi_m = project_from_slave_to_master(master_element, X_s, n_s, time)
                 N2 = vec(get_basis(master_element, xi_m, time))
-                X_m = N2*X2 
+                X_m = N2*X2
                 De += w*Phi*N1'
                 Me += w*Phi*N2'
                 if props.adjust
@@ -194,7 +194,7 @@ function assemble!(problem::Problem{Mortar}, time::Float64, ::Type{Val{1}}, ::Ty
             # add contribution to contact virtual work
             sdofs = get_gdofs(problem, slave_element)
             mdofs = get_gdofs(problem, master_element)
-            
+
             for i=1:field_dim
                 lsdofs = sdofs[i:field_dim:end]
                 lmdofs = mdofs[i:field_dim:end]
@@ -210,4 +210,3 @@ function assemble!(problem::Problem{Mortar}, time::Float64, ::Type{Val{1}}, ::Ty
     end # slave elements done, contact virtual work ready
 
 end
-
diff --git a/src/solvers.jl b/src/solvers.jl
index ba4c4dc..5e8b66c 100644
--- a/src/solvers.jl
+++ b/src/solvers.jl
@@ -10,19 +10,21 @@ type Solver{S<:AbstractSolver}
     norms :: Vector{Tuple}       # solution norms for convergence studies
     ndofs :: Int                 # number of degrees of freedom in problem
     xdmf :: Nullable{Xdmf}       # input/output handle
+    initialized :: Bool
+    u :: Vector{Float64}
+    la :: Vector{Float64}
     properties :: S
 end
 
 
 function Solver{S<:AbstractSolver}(::Type{S}, name="solver", properties...)
     variant = S(properties...)
-    solver = Solver{S}(name, 0.0, [], [], 0, nothing, variant)
+    solver = Solver{S}(name, 0.0, [], [], 0, nothing, false, [], [], variant)
     return solver
 end
 
 function Solver{S<:AbstractSolver}(::Type{S}, problems::Problem...)
-    variant = S()
-    solver = Solver{S}("$(S)Solver", 0.0, [], [], 0, nothing, variant)
+    solver = Solver(S, "$(S)Solver")
     push!(solver.problems, problems...)
     return solver
 end
@@ -237,13 +239,13 @@ Solve linear system using LDLt factorization (SuiteSparse). This version
 requires that final system is symmetric and positive definite, so boundary
 conditions are first eliminated before solution.
 """
-function solve!(K, C1, C2, D, f, g, u, la, ::Type{Val{1}}; F=nothing, debug=false)
+function solve!(solver::Solver, K, C1, C2, D, f, g, u, la, ::Type{Val{1}}; debug=false)
 
-    nnz(D) == 0 || return F, false
+    nnz(D) == 0 || return false
     nz = get_nonzero_rows(C2)
     B = get_nonzero_rows(C2')
     # C2^-1 exists or this doesn't work
-    length(nz) == length(B) || return F, false
+    length(nz) == length(B) || return false
 
     A = get_nonzero_rows(K)
     I = setdiff(A, B)
@@ -270,38 +272,30 @@ function solve!(K, C1, C2, D, f, g, u, la, ::Type{Val{1}}; F=nothing, debug=fals
     end
 
     # solve interior domain using LDLt factorization
-    if F == nothing
-        F = ldltfact(K[I,I])
-    end
+    F = ldltfact(K[I,I])
     u[I] = F \ (f[I] - K[I,B]*u[B])
     # solve lambda
     la[B] = lufact(C1[B,nz]) \ full(f[B] - K[B,I]*u[I] - K[B,B]*u[B])
 
-    return F, true
+    return true
 end
 
 """
 Solve linear system using LU factorization (UMFPACK). This version solves
 directly the saddle point problem without elimination of boundary conditions.
 """
-function solve!(K, C1, C2, D, f, g, u, la, ::Type{Val{2}}; F=nothing)
+function solve!(solver::Solver, K, C1, C2, D, f, g, u, la, ::Type{Val{2}})
     # construct global system Ax = b and solve using lufact (UMFPACK)
     A = [K C1'; C2  D]
     b = [f; g]
-    nz = get_nonzero_rows(A)
-    x = zeros(length(b))
-    if F == nothing
-        F = lufact(A[nz,nz])
-    end
-    x[nz] = F \ full(b[nz])
-    ndofs = size(K, 1)
-    u[:] = x[1:ndofs]
-    la[:] = x[ndofs+1:end]
-    return F, true
+    x = lufact(A) \ full(b)
+    u[:] = x[1:solver.ndofs]
+    la[:] = x[solver.ndofs+1:end]
+    return true
 end
 
 """ Default linear system solver for solver. """
-function solve_linear_system(solver::Solver; F=nothing, empty_assemblies_before_solution=true, show_info=true)
+function solve!(solver::Solver; empty_assemblies_before_solution=true, show_info=true)
     show_info && info("Solving problems ...")
     t0 = Base.time()
 
@@ -317,6 +311,11 @@ function solve_linear_system(solver::Solver; F=nothing, empty_assemblies_before_
     K = 1/2*(K + K')
     M = 1/2*(M + M')
 
+    nz = ones(solver.ndofs)
+    nz[get_nonzero_rows(C2)] = 0.0
+    nz[get_nonzero_rows(D)] = 0.0
+    D += spdiagm(nz)
+
     # free up some memory before solution
     for problem in get_problems(solver)
         if empty_assemblies_before_solution
@@ -327,22 +326,25 @@ function solve_linear_system(solver::Solver; F=nothing, empty_assemblies_before_
         gc()
     end
 
-    u = zeros(solver.ndofs)
-    la = zeros(solver.ndofs)
-
+    ndofs = solver.ndofs
+    u = zeros(ndofs)
+    la = zeros(ndofs)
     status = false
     i = 0
     for i in [1, 2]
-        F, status = solve!(K, C1, C2, D, f, g, u, la, Val{i}; F=F)
+        status = solve!(solver, K, C1, C2, D, f, g, u, la, Val{i})
         status && break
     end
     status || error("Failed to solve linear system!")
-
     t1 = round(Base.time()-t0, 2)
     norms = (norm(u), norm(la))
     show_info && info("Solved problems in $t1 seconds using solver $i. Solution norms = $norms.")
     push!(solver.norms, norms)
-    return F, u, la
+
+    solver.u = u
+    solver.la = la
+
+    return
 end
 
 """ Default assembler for solver. """
@@ -398,7 +400,11 @@ end
 
 """ Default initializer for solver. """
 function initialize!(solver::Solver; show_info=true)
-    show_info && info("Initializing problems ...")
+    if solver.initialized
+        show_info && info("initialize!(): solver already initialized")
+        return
+    end
+    show_info && info("Initializing solver ...")
     problems = get_problems(solver)
     length(problems) != 0 || error("Empty solver, add problems to solver using push!")
     t0 = Base.time()
@@ -419,16 +425,18 @@ function initialize!(solver::Solver; show_info=true)
     info("Total number of nodes in problems: $nnodes")
     maxdof = maximum(nnodes)*field_dim
     info("# of max dof (=size of solution vector) is $maxdof")
-    u = zeros(maxdof)
-    la = zeros(maxdof)
+    solver.u = zeros(maxdof)
+    solver.la = zeros(maxdof)
     # TODO: this could be used to initialize elements too...
+    # TODO: cannot initialize to zero always, construct vector from elements.
     for problem in problems
-        problem.assembly.u = u
-        problem.assembly.la = la
+        problem.assembly.u = zeros(maxdof)
+        problem.assembly.la = zeros(maxdof)
         # initialize(problem, ....)
     end
     t1 = round(Base.time()-t0, 2)
-    show_info && info("Initialized problems in $t1 seconds.")
+    show_info && info("Initialized solver in $t1 seconds.")
+    solver.initialized = true
 end
 
 function get_all_elements(solver::Solver)
@@ -472,7 +480,10 @@ function get_temporal_collection(xdmf::Xdmf)
 end
 
 """ Default update for solver. """
-function update!(solver::Solver, u::Vector, la::Vector; show_info=true)
+function update!{S}(solver::Solver{S}; show_info=true)
+    u = solver.u
+    la = solver.la
+
     show_info && info("Updating problems ...")
     t0 = Base.time()
 
@@ -487,88 +498,105 @@ function update!(solver::Solver, u::Vector, la::Vector; show_info=true)
 
     # if io is attached to solver, update hdf / xml also
     if !isnull(solver.xdmf)
-        xdmf = get(solver.xdmf)
-        temporal_collection = get_temporal_collection(xdmf)
-        frame = new_child(temporal_collection, "Grid")
-        new_child(frame, "Time", Dict("Value" => solver.time))
-
-        # save geometry
-        X = solver("geometry", solver.time)
-        node_ids = sort(collect(keys(X)))
-        geometry = hcat([X[nid] for nid in node_ids]...)
-        ndim, nnodes = size(geometry)
-        geom_type = ndim == 2 ? "XY" : "XYZ"
-        dataitem = new_dataitem(xdmf, "/Node IDs", node_ids)
-        geom = new_child(frame, "Geometry", Dict("Type" => geom_type))
-        dataitem = new_dataitem(xdmf, "/Geometry", geometry)
-        add_child(geom, dataitem)
-
-        # save topology
-        all_elements = get_all_elements(solver)
-        nelements = length(all_elements)
-        element_types = unique(map(get_element_type, all_elements))
-
-        xdmf_element_mapping = Dict(
-            "Seg2" => "Polyline",
-            "Tri3" => "Triangle",
-            "Quad4" => "Quadrilateral",
-            "Tet4" => "Tetrahedron",
-            "Pyramid5" => "Pyramid",
-            "Wedge6" => "Wedge",
-            "Hex8" => "Hexahedron",
-            "Seg3" => "Edge_3",
-            "Tri6" => "Tri_6",
-            "Quad8" => "Quad_8",
-            "Tet10" => "Tet_10",
-            "Pyramid13" => "Pyramid_13",
-            "Wedge15" => "Wedge_15",
-            "Hex20" => "Hex_20")
-
-        for element_type in element_types
-            elements = filter_by_element_type(element_type, all_elements)
-            sort!(elements, by=get_element_id)
-            element_ids = map(get_element_id, elements)
-            element_conn = map(get_connectivity, elements)
-            element_conn = transpose(hcat(element_conn...)) - 1
-            element_code = split(string(element_type), ".")[end]
-            dataitem = new_dataitem(xdmf, "/Topology/$element_code/Element IDs", element_ids)
-            dataitem = new_dataitem(xdmf, "/Topology/$element_code/Connectivity", element_conn)
-            topology = new_child(frame, "Topology")
-            set_attribute(topology, "TopologyType", xdmf_element_mapping[element_code])
-            set_attribute(topology, "NumberOfElements", length(elements))
-            add_child(topology, dataitem)
-        end
-
-        # save solved fields
-        time = "Time $(solver.time)"
-        unknown_field_name = get_unknown_field_name(solver)
-        U = solver(unknown_field_name, solver.time)
-        node_ids2 = sort(collect(keys(U)))
-
-        @assert node_ids == node_ids2
-        ndim = length(U[first(node_ids)])
-        field_type = ndim == 1 ? "Scalar" : "Vector"
-        field_center = "Node"
-        if ndim == 2
-            for nid in node_ids
-                U[nid] = [U[nid]; 0.0]
-            end
-            ndim = 3
-        end
-        U = hcat([U[nid] for nid in node_ids]...)
-        unknown_field_name = ucfirst(unknown_field_name)
-        dataitem = new_dataitem(xdmf, "/Results/$time/Nodal Fields/$unknown_field_name", U)
-        attribute = new_child(frame, "Attribute")
-        set_attribute(attribute, "Name", unknown_field_name)
-        set_attribute(attribute, "Center", field_center)
-        add_child(attribute, dataitem)
-        save!(xdmf)
+        update_xdmf!(solver)
     end
 
     t1 = round(Base.time()-t0, 2)
     show_info && info("Updated problems in $t1 seconds.")
 end
 
+function update_xdmf!{S}(solver::Solver{S}; show_info=true)
+    xdmf = get(solver.xdmf)
+    temporal_collection = get_temporal_collection(xdmf)
+
+    frame = new_element("Grid")
+    new_child(frame, "Time", Dict("Value" => solver.time))
+
+    # save geometry
+    X = solver("geometry", solver.time)
+    node_ids = sort(collect(keys(X)))
+    geometry = hcat([X[nid] for nid in node_ids]...)
+    ndim, nnodes = size(geometry)
+    geom_type = ndim == 2 ? "XY" : "XYZ"
+    dataitem = new_dataitem(xdmf, "/Node IDs", node_ids)
+    geom = new_child(frame, "Geometry", Dict("Type" => geom_type))
+    dataitem = new_dataitem(xdmf, "/Geometry", geometry)
+    add_child(geom, dataitem)
+
+    # save topology
+    all_elements = get_all_elements(solver)
+    nelements = length(all_elements)
+    element_types = unique(map(get_element_type, all_elements))
+
+    xdmf_element_mapping = Dict(
+        "Seg2" => "Polyline",
+        "Tri3" => "Triangle",
+        "Quad4" => "Quadrilateral",
+        "Tet4" => "Tetrahedron",
+        "Pyramid5" => "Pyramid",
+        "Wedge6" => "Wedge",
+        "Hex8" => "Hexahedron",
+        "Seg3" => "Edge_3",
+        "Tri6" => "Tri_6",
+        "Quad8" => "Quad_8",
+        "Tet10" => "Tet_10",
+        "Pyramid13" => "Pyramid_13",
+        "Wedge15" => "Wedge_15",
+        "Hex20" => "Hex_20")
+
+    for element_type in element_types
+        elements = filter_by_element_type(element_type, all_elements)
+        sort!(elements, by=get_element_id)
+        element_ids = map(get_element_id, elements)
+        element_conn = map(get_connectivity, elements)
+        element_conn = transpose(hcat(element_conn...)) - 1
+        element_code = split(string(element_type), ".")[end]
+        dataitem = new_dataitem(xdmf, "/Topology/$element_code/Element IDs", element_ids)
+        dataitem = new_dataitem(xdmf, "/Topology/$element_code/Connectivity", element_conn)
+        topology = new_child(frame, "Topology")
+        set_attribute(topology, "TopologyType", xdmf_element_mapping[element_code])
+        set_attribute(topology, "NumberOfElements", length(elements))
+        add_child(topology, dataitem)
+    end
+
+    # save solved fields
+    unknown_field_name = get_unknown_field_name(solver)
+    U = solver(unknown_field_name, solver.time)
+    node_ids2 = sort(collect(keys(U)))
+
+    @assert node_ids == node_ids2
+    ndim = length(U[first(node_ids)])
+    field_type = ndim == 1 ? "Scalar" : "Vector"
+    field_center = "Node"
+    if ndim == 2
+        for nid in node_ids
+            U[nid] = [U[nid]; 0.0]
+        end
+        ndim = 3
+    end
+    U = hcat([U[nid] for nid in node_ids]...)
+    unknown_field_name = ucfirst(unknown_field_name)
+    time = solver.time
+    path = ""
+    if S == Nonlinear
+        iteration = solver.properties.iteration
+        path = "/Results/Time $time/Iteration $iteration/Nodal Fields/$unknown_field_name"
+    elseif S == Linear
+        path = "/Results/Time $time/Nodal Fields/$unknown_field_name"
+    end
+    dataitem = new_dataitem(xdmf, path, U)
+    attribute = new_child(frame, "Attribute")
+    set_attribute(attribute, "Name", unknown_field_name)
+    set_attribute(attribute, "Center", field_center)
+    set_attribute(attribute, "AttributeType", field_type)
+    add_child(attribute, dataitem)
+    if (S == Linear) || ((S == Nonlinear) && has_converged(solver))
+        add_child(temporal_collection, frame)
+    end
+    save!(xdmf)
+end
+
+
 ### Nonlinear quasistatic solver
 
 type Nonlinear <: AbstractSolver
@@ -580,7 +608,7 @@ type Nonlinear <: AbstractSolver
 end
 
 function Nonlinear()
-    solver = Nonlinear(0, 1, 20, 5.0e-5, true)
+    solver = Nonlinear(0, 1, 10, 5.0e-5, true)
     return solver
 end
 
@@ -646,12 +674,10 @@ function call(solver::Solver{Nonlinear}; show_info=true)
 
         # 2.1 update linearized assemblies
         assemble!(solver)
-
         # 2.2 call solver for linearized system
-        F, u, la = solve_linear_system(solver)
-
+        solve!(solver)
         # 2.3 update solution back to elements
-        update!(solver, u, la)
+        update!(solver)
 
         # 2.4 check convergence
         if has_converged(solver)
@@ -713,7 +739,7 @@ function assemble!(solver::Solver{Linear}; show_info=true)
     show_info && info("Assembled $nproblems problems in $t1 seconds. ndofs = $ndofs.")
 end
 
-function call(solver::Solver{Linear}; F=nothing, show_info=true, return_factorization=true)
+function call(solver::Solver{Linear}; show_info=true)
     t0 = Base.time()
     show_info && info(repeat("-", 80))
     show_info && info("Starting linear solver")
@@ -721,13 +747,10 @@ function call(solver::Solver{Linear}; F=nothing, show_info=true, return_factoriz
     show_info && info(repeat("-", 80))
     initialize!(solver)
     assemble!(solver)
-    F, u, la = solve_linear_system(solver; F=F, empty_assemblies_before_solution=false)
-    update!(solver, u, la)
+    solve!(solver)
+    update!(solver)
     t1 = round(Base.time()-t0, 2)
     show_info && info("Linear solver ready in $t1 seconds.")
-    if return_factorization
-        return F
-    end
 end
 
 """ Convenience function to call linear solver. """
diff --git a/test/test_abaqus_verification.jl b/test/test_abaqus_verification.jl
index 5225af1..9316084 100644
--- a/test/test_abaqus_verification.jl
+++ b/test/test_abaqus_verification.jl
@@ -6,6 +6,7 @@ using JuliaFEM.Preprocess
 using JuliaFEM.Postprocess
 using JuliaFEM.Abaqus
 using JuliaFEM.Testing
+using DataFrames
 
 # to turn on automatic file download, set
 # ENV["ABAQUS_DOWNLOAD_URL"] = "http://<domain>:2080/v2016/books/eif"
@@ -47,8 +48,32 @@ end
       @testset "1.3.3 Three-dimensional solid elements" begin
         @testset "C3D8 elements." begin
           abaqus_run_test("ec38sfs2") || return
+          res = abaqus_open_results("ec38sfs2")
+
+node_output1 = wsv"""
+NODE           U1          U2          U3          COOR1       COOR2       COOR3
+1     -2.0000E-33 -2.0000E-33 -2.0000E-33          0.000       0.000       0.000
+2     -2.6667E-05 -1.0000E-33 -1.7333E-04          2.000       0.000       0.000
+3     -2.0000E-04 -2.6667E-05 -1.7333E-04          2.000       2.000       0.000
+4     -1.7333E-04 -2.6667E-05 -1.0000E-33          0.000       2.000       0.000
+5     -3.6777E-48 -8.6667E-05 -1.3333E-05          0.000       0.000       1.000
+6     -2.6667E-05 -8.6667E-05 -1.8667E-04          2.000       0.000       1.000
+7     -2.0000E-04 -1.1333E-04 -1.8667E-04          2.000       2.000       1.000
+8     -1.7333E-04 -1.1333E-04 -1.3333E-05          0.000       2.000       1.000
+"""
+
+node_output_2 = wsv"""
+NODE          RF1         RF2         RF3            CF1         CF2         CF3
+1        1500.000    1500.000    1000.000          0.000       0.000       0.000
+2           0.000     500.000       0.000       1500.000       0.000       0.000
+3           0.000       0.000       0.000        500.000     500.000   -1000.000
+4           0.000       0.000       0.000        500.000    1500.000       0.000
+5        -500.000       0.000       0.000          0.000    -500.000    1000.000
+6           0.000       0.000       0.000       -500.000   -1500.000       0.000
+7           0.000       0.000       0.000      -1500.000   -1500.000   -1000.000
+8           0.000       0.000       0.000      -1500.000    -500.000       0.000
+"""
           #= to check also results:
-          xdmf = abaqus_open_results("ec38sfs2")
           side, opts = read_result(xdmf, "SECTION/side")
           @test isapprox(side["SOFM"], 3464.0)
           @test isapprox(side["SOF1"], 2000.0)
@@ -68,4 +93,3 @@ end
     end
   end
 end
-
diff --git a/test/test_elasticity_2d_linear_with_surface_load.jl b/test/test_elasticity_2d_linear_with_surface_load.jl
index 6271663..ce2cf98 100644
--- a/test/test_elasticity_2d_linear_with_surface_load.jl
+++ b/test/test_elasticity_2d_linear_with_surface_load.jl
@@ -32,7 +32,7 @@ using JuliaFEM.Testing
     update!(block.elements, "displacement load 2", 576.0)
 
     # traction
-    traction = Problem(Elasticity, "BLOCK", 2)
+    traction = Problem(Elasticity, "TRACTION", 2)
     traction.properties.formulation = :plane_stress
     traction.properties.finite_strain = false
     traction.properties.geometric_stiffness = false
@@ -48,8 +48,6 @@ using JuliaFEM.Testing
     update!(bc_sym_13, "displacement 2", 0.0)
 
     solver = LinearSolver(block, traction, bc_sym_23, bc_sym_13)
-#   assemble!(solver)
-#   dump(full(bc_sym_23.assembly.C1))
     solver()
 
     info("u = ", block.assembly.u)
@@ -98,6 +96,19 @@ using JuliaFEM.Testing
     S = solver(DataFrame, 0.0, Val{:S})
     println(S)
 
+    info(solver("displacement", 0.0))
+    solver()
+    info(solver("displacement", 0.0))
+    u = solver("displacement", 0.0)[3]
+    info("u3 = $u")
+    @test isapprox(u, u3_expected)
+
+    info("calling nonlinear solver")
+    solver2 = NonlinearSolver(block, traction, bc_sym_23, bc_sym_13)
+    solver2()
+    u = solver2("displacement", 0.0)[3]
+    info("nlsolver u3 = $u, expected = $u3_expected")
+    @test isapprox(u, u3_expected; rtol=1.0e-5)
 end
 
 #= TODO: to other file
diff --git a/test/test_elements_add_fields.jl b/test/test_elements_add_fields.jl
index a9c4b81..b85db1a 100644
--- a/test/test_elements_add_fields.jl
+++ b/test/test_elements_add_fields.jl
@@ -4,4 +4,15 @@
 using JuliaFEM
 using JuliaFEM.Testing
 
-
+@testset "dict field" begin
+    el = Element(Seg2, 1, [1, 2])
+    X = Dict{Int64, Vector{Float64}}(1 => [0.0, 0.0], 2 => [1.0, 0.0], 3 => [0.5, 0.5])
+    f = Field(X)
+    debug("field = $f")
+    #update!(el, "geometry", X)
+    el["geometry"] = f
+    @test isapprox(el("geometry")[1], [0.0, 0.0])
+    @test isapprox(el("geometry", 0.0)[1], [0.0, 0.0])
+    @test isapprox(el("geometry", 0.0)[3], [0.5, 0.5])
+    @test isapprox(el("geometry", [0.0], 0.0), [0.5, 0.0])
+end
diff --git a/test/test_fields.jl b/test/test_fields.jl
index 8426ea8..667a616 100644
--- a/test/test_fields.jl
+++ b/test/test_fields.jl
@@ -3,8 +3,12 @@
 
 using JuliaFEM
 using JuliaFEM.Testing
+using Logging
+Logging.configure(level=DEBUG)
 
-@testset "test updating time dependent fields" begin
+@testset "create and manipulate fields" begin
+
+@testset "updating time dependent fields" begin
     f = Field(0.0 => 1.0)
     @test last(f).time == 0.0
     @test last(f).data == 1.0
@@ -18,16 +22,41 @@ using JuliaFEM.Testing
     @test length(f) == 2
 end
 
-@testset "test updating time invariant fields" begin
+@testset "updating time invariant fields" begin
     f = Field(1.0)
     @test f.data == 1.0
     update!(f, 2.0)
     @test f.data == 2.0
 end
 
-@testset "test field defined using function" begin
+@testset "field defined using function" begin
     g(xi, t) = xi[1]*t
     f = Field(g)
     v = f([1.0], 2.0)
     @test isapprox(v, 2.0)
 end
+
+@testset "dictionary fields" begin
+    f1 = Dict{Int64, Vector{Float64}}(1 => [0.0, 0.0], 2 => [0.0, 0.0])
+    f2 = Dict{Int64, Vector{Float64}}(1 => [1.0, 1.0], 2 => [1.0, 1.0])
+    f = Field(0.0 => f1, 1.0 => f2)
+    debug("field = $f")
+    @test isa(f, DVTV)
+    @test isapprox(f(0.0)[1], [0.0, 0.0])
+    @test isapprox(f(1.0)[2], [1.0, 1.0])
+
+    f = Field(0.0 => f1)
+    update!(f, 1.0 => f2)
+    @test isa(f, DVTV)
+    @test isapprox(f(0.0)[1], [0.0, 0.0])
+    @test isapprox(f(1.0)[2], [1.0, 1.0])
+
+    f = Field(f1)
+    @test isapprox(f(0.0)[1], [0.0, 0.0])
+    @test isapprox(f[1], [0.0, 0.0])
+
+    f = Field(f1)
+    @test isa(f, DVTI)
+end
+
+end
diff --git a/test/test_heat_4.jl b/test/test_heat_4.jl
new file mode 100644
index 0000000..2267abf
--- /dev/null
+++ b/test/test_heat_4.jl
@@ -0,0 +1,92 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
+
+using JuliaFEM
+using JuliaFEM.Testing
+using DataFrames
+
+@testset "two increments, linear solver" begin
+    X = Dict{Int, Vector{Float64}}(
+        1 => [0.0,0.0],
+        2 => [1.0,0.0],
+        3 => [1.0,1.0],
+        4 => [0.0,1.0])
+    element = Element(Quad4, [1, 2, 3, 4])
+    update!(element, "geometry", X)
+    update!(element, "temperature thermal conductivity", 6.0)
+    update!(element, "temperature load", 0.0 => 12.0)
+    update!(element, "temperature load", 1.0 => 24.0)
+    problem = Problem(Heat, "one element heat problem", 1)
+    problem.properties.formulation = "2D"
+    push!(problem, element)
+    boundary_element = Element(Seg2, [1, 2])
+    update!(boundary_element, "geometry", X)
+    update!(boundary_element, "temperature 1", 0.0)
+    bc = Problem(Dirichlet, "fixed", 1, "temperature")
+    push!(bc, boundary_element)
+    solver = Solver(Linear, problem, bc)
+
+    solver.time = 0.0
+
+    empty!(problem.assembly)
+    solver()
+    @test isapprox(solver("temperature", 0.0)[3], 1.0)
+
+    empty!(problem.assembly)
+    solver()
+    @test isapprox(solver("temperature", 0.0)[3], 1.0)
+
+    solver.time = 1.0
+
+    empty!(problem.assembly)
+    solver()
+    @test isapprox(solver("temperature", 1.0)[3], 2.0)
+
+    empty!(problem.assembly)
+    solver()
+    @test isapprox(solver("temperature", 1.0)[3], 2.0)
+
+end
+
+@testset "two increments, nonlinear solver" begin
+    X = Dict{Int, Vector{Float64}}(
+        1 => [0.0,0.0],
+        2 => [1.0,0.0],
+        3 => [1.0,1.0],
+        4 => [0.0,1.0])
+    element = Element(Quad4, [1, 2, 3, 4])
+    update!(element, "geometry", X)
+    update!(element, "temperature thermal conductivity", 6.0)
+    update!(element, "temperature load", 0.0 => 12.0)
+    update!(element, "temperature load", 1.0 => 24.0)
+    problem = Problem(Heat, "one element heat problem", 1)
+    problem.properties.formulation = "2D"
+    push!(problem, element)
+    boundary_element = Element(Seg2, [1, 2])
+    update!(boundary_element, "geometry", X)
+    update!(boundary_element, "temperature 1", 0.0)
+    bc = Problem(Dirichlet, "fixed", 1, "temperature")
+    push!(bc, boundary_element)
+    solver = Solver(Nonlinear, problem, bc)
+
+    solver.time = 0.0
+
+    empty!(problem.assembly)
+    solver()
+    @test isapprox(solver("temperature", 0.0)[3], 1.0)
+
+    empty!(problem.assembly)
+    solver()
+    @test isapprox(solver("temperature", 0.0)[3], 1.0)
+
+    solver.time = 1.0
+
+    empty!(problem.assembly)
+    solver()
+    @test isapprox(solver("temperature", 1.0)[3], 2.0)
+
+    empty!(problem.assembly)
+    solver()
+    @test isapprox(solver("temperature", 1.0)[3], 2.0)
+
+end
diff --git a/test/test_io.jl b/test/test_io.jl
index 077519d..a997cac 100644
--- a/test/test_io.jl
+++ b/test/test_io.jl
@@ -55,7 +55,52 @@ end
     @test isapprox(read(xdmf, "/Domain/Grid/Grid[2]/Geometry/DataItem"), [1.0, 2.0])
 end
 
-@testset "save results to disk" begin
+@testset "higher level xdmf" begin
+    e1 = Element(Quad4, 1, [1, 2, 3, 4])
+    e2 = Element(Quad4, 2, [5, 6, 7, 8])
+    p1 = Problem(Elasticity, "Body 1", 2)
+    p2 = Problem(Elasticity, "Body 2", 2)
+    push!(p1, e1)
+    push!(p2, e2)
+    #update!(p1)
+
+    e3 = Element(Seg2, 3, [1, 2])
+    e4 = Element(Seg2, 4, [3, 4])
+    e5 = Element(Seg2, 5, [5, 6])
+    p3 = Problem(Dirichlet, "Fixed BC", 2, "displacement")
+    p4 = Problem(Contact, "Contact between bodies 1 and 2", 2, "displacement")
+    push!(p3, e3)
+    push!(p4, e4)
+    X = Dict{Int64, Vector{Float64}}(
+        1 => [0.0, 0.0],
+        2 => [1.0, 0.0],
+        3 => [1.0, 1.0],
+        4 => [0.0, 1.0],
+        5 => [0.0, 2.0],
+        6 => [1.0, 2.0],
+        7 => [1.0, 3.0],
+        8 => [0.0, 3.0])
+    u = Dict{Int64, Vector{Float64}}(
+        1 => [0.1, 0.1],
+        2 => [0.1, 0.1],
+        3 => [0.1, 0.1],
+        4 => [0.1, 0.1],
+        5 => [0.1, 0.1],
+        6 => [0.1, 0.1],
+        7 => [0.1, 0.1],
+        8 => [0.1, 0.1])
+    n = Dict{Int64, Vector{Float64}}(
+        3 => [0.0, 1.0],
+        4 => [0.0, 1.0])
+    R = Dict{Int64, Vector{Float64}}(
+        1 => [0.0, 1.0],
+        2 => [0.0, 1.0])
+    update!(e4, "master elements", [e3])
+end
+
+#=
+
+@testset "save results to disk, linear solver" begin
     X = Dict{Int, Vector{Float64}}(
         1 => [0.0,0.0],
         2 => [1.0,0.0],
@@ -65,7 +110,7 @@ end
     update!(element, "geometry", X)
     update!(element, "temperature thermal conductivity", 6.0)
     update!(element, "temperature load", 0.0 => 12.0)
-    update!(element, "temperature load", 1.0 => 18.0)
+    update!(element, "temperature load", 1.0 => 24.0)
     problem = Problem(Heat, "one element heat problem", 1)
     problem.properties.formulation = "2D"
     push!(problem, element)
@@ -74,8 +119,9 @@ end
     update!(boundary_element, "temperature 1", 0.0)
     bc = Problem(Dirichlet, "fixed", 1, "temperature")
     push!(bc, boundary_element)
+    xdmf = Xdmf()
     solver = Solver(Linear, problem, bc)
-    solver.xdmf = Xdmf()
+    solver.xdmf = xdmf
 
     solver.time = 0.0
     solver()
@@ -88,7 +134,6 @@ end
     info(element("temperature load", [0.0, 0.0], 0.0))
     info(element("temperature load", [0.0, 0.0], 1.0))
 
-    xdmf = get(solver.xdmf)
     info("h5 file = $(h5file(xdmf))")
     E = read(xdmf.hdf, "/Topology/Quad4/Element IDs")
     C = read(xdmf.hdf, "/Topology/Quad4/Connectivity")
@@ -101,7 +146,7 @@ end
     @test isapprox(N, [1, 2, 3, 4])
     X_expected = [0.0 0.0; 1.0 0.0; 1.0 1.0; 0.0 1.0]'
     T1_expected = [0.0 0.0 1.0 1.0]
-    T2_expected = [0.0 0.0 0.5 0.5]
+    T2_expected = [0.0 0.0 2.0 2.0]
     @test isapprox(X, X_expected)
     @test isapprox(T1, T1_expected)
     @test isapprox(T2, T2_expected)
@@ -115,3 +160,62 @@ end
     @test isapprox(read(xdmf, "/Domain/Grid/Grid[end]/Time/Value"), 1.0)
     @test isapprox(read(xdmf, "/Domain/Grid/Grid[end]/Topology/DataItem"), [0 1 2 3])
 end
+
+@testset "save results to disk, nonlinear solver" begin
+    X = Dict{Int, Vector{Float64}}(
+        1 => [0.0,0.0],
+        2 => [1.0,0.0],
+        3 => [1.0,1.0],
+        4 => [0.0,1.0])
+    element = Element(Quad4, [1, 2, 3, 4])
+    update!(element, "geometry", X)
+    update!(element, "temperature thermal conductivity", 6.0)
+    update!(element, "temperature load", 0.0 => 12.0)
+    update!(element, "temperature load", 1.0 => 24.0)
+    problem = Problem(Heat, "one element heat problem", 1)
+    problem.properties.formulation = "2D"
+    push!(problem, element)
+    boundary_element = Element(Seg2, [1, 2])
+    update!(boundary_element, "geometry", X)
+    update!(boundary_element, "temperature 1", 0.0)
+    bc = Problem(Dirichlet, "fixed", 1, "temperature")
+    push!(bc, boundary_element)
+    solver = Solver(Nonlinear, problem, bc)
+    solver.xdmf = Xdmf()
+
+    solver.time = 0.0
+    solver()
+    solver.time = 1.0
+    solver()
+
+    xdmf = get(solver.xdmf)
+    info("h5 file = $(h5file(xdmf))")
+    E = read(xdmf.hdf, "/Topology/Quad4/Element IDs")
+    C = read(xdmf.hdf, "/Topology/Quad4/Connectivity")
+    N = read(xdmf.hdf, "/Node IDs")
+    X = read(xdmf.hdf, "/Geometry")
+    T11 = read(xdmf.hdf, "/Results/Time 0.0/Iteration 1/Nodal Fields/Temperature")
+    T12 = read(xdmf.hdf, "/Results/Time 0.0/Iteration 2/Nodal Fields/Temperature")
+    T21 = read(xdmf.hdf, "/Results/Time 1.0/Iteration 1/Nodal Fields/Temperature")
+    T22 = read(xdmf.hdf, "/Results/Time 1.0/Iteration 2/Nodal Fields/Temperature")
+    X_expected = [
+        0.0 0.0
+        1.0 0.0
+        1.0 1.0
+        0.0 1.0]
+    T1_expected = [0.0 0.0 1.0 1.0]
+    T2_expected = [0.0 0.0 2.0 2.0]
+    @test isapprox(T12, T1_expected)
+    @test isapprox(T22, T2_expected)
+    @test isapprox(read(xdmf, "/Domain/Grid/Grid/Time/Value"), 0.0)
+    @test read(xdmf, "/Domain/Grid/Grid/Geometry/Type") == "XY"
+    @test isapprox(read(xdmf, "/Domain/Grid/Grid/Geometry/DataItem"), X_expected')
+    @test isapprox(read(xdmf, "/Domain/Grid/Grid/Topology/DataItem"), [0 1 2 3])
+    @test isapprox(read(xdmf, "/Domain/Grid/Grid/Topology[@TopologyType=Polyline]/DataItem"), [0 1])
+    @test isapprox(read(xdmf, "/Domain/Grid/Grid[1]/Attribute[@Name=Temperature]/DataItem"), T1_expected)
+    @test isapprox(read(xdmf, "/Domain/Grid/Grid[2]/Attribute[@Name=Temperature]/DataItem"), T2_expected)
+    @test isapprox(read(xdmf, "/Domain/Grid/Grid[end]/Time/Value"), 1.0)
+    @test isapprox(read(xdmf, "/Domain/Grid/Grid[end]/Topology/DataItem"), [0 1 2 3])
+end
+
+=#
diff --git a/test/test_mortar_2d_mesh_tie.jl b/test/test_mortar_2d_mesh_tie.jl
index fc38a8c..64f7d21 100644
--- a/test/test_mortar_2d_mesh_tie.jl
+++ b/test/test_mortar_2d_mesh_tie.jl
@@ -88,13 +88,13 @@ end
     push!(solver, upper, lower, bc_upper, bc_lower, interface)
     solver()
 
-    interface_norm = norm(interface.assembly)
+    #interface_norm = norm(interface.assembly)
     # for bi-orthogonal:
     #interface_norm_expected = [0.0, 0.0, 0.0, 0.0, 0.0, 0.44870723441585775, 0.44870723441585775, 0.0, 0.0, 0.0]
-    interface_norm_expected = [0.0, 0.0, 0.0, 0.0, 0.0, 0.39361633468943247, 0.39361633468943247, 0.0, 0.0, 0.0]
-    info("Interface norm:          $interface_norm")
-    info("Interface norm expected: $interface_norm_expected")
-    @test isapprox(interface_norm, interface_norm_expected)
+    #interface_norm_expected = [0.0, 0.0, 0.0, 0.0, 0.0, 0.39361633468943247, 0.39361633468943247, 0.0, 0.0, 0.0]
+    #info("Interface norm:          $interface_norm")
+    #info("Interface norm expected: $interface_norm_expected")
+    #@test isapprox(interface_norm, interface_norm_expected)
 
     T_upper = first(bc_upper.elements)("temperature", [0.0], 0.0)
     T_lower = first(bc_lower.elements)("temperature", [0.0], 0.0)
@@ -285,4 +285,3 @@ end
     interface = solver["interface between upper and lower block"]
     @test isapprox(norm(interface.assembly.u), 0.34318800698017704)
 end
-
diff --git a/test/test_problem.jl b/test/test_problem.jl
index 035d425..aef7c74 100644
--- a/test/test_problem.jl
+++ b/test/test_problem.jl
@@ -13,7 +13,7 @@ using JuliaFEM.Testing
     # one timestep in field "temperature"
     @test length(el["temperature"]) == 1
     # this way we access to field at default time t=0.0, it's different than ^!
-    @test length(el("temperature")) == 2 
+    @test length(el("temperature")) == 2
     # length of single increment
     @test length(el("temperature", 0.0)) == 2
     @test length(last(el, "temperature").data) == 2
@@ -27,7 +27,7 @@ end
     @test haskey(el, "displacement")
     @test length(el["displacement"]) == 1
     # this way we access to field at default time t=0.0, it's different than ^!
-    @test length(el("displacement")) == 2 
+    @test length(el("displacement")) == 2
     # length of single increment
     @test length(el("displacement", 0.0)) == 2
     @test length(last(el, "displacement").data) == 2
@@ -41,3 +41,45 @@ end
     @test haskey(el, "reaction force")
     @test haskey(el, "temperature")
 end
+
+#=
+@testset "dict field depending from problems" begin
+    p1 = Problem(Elasticity, "Body 1", 2)
+    p2 = Problem(Elasticity, "Body 2", 2)
+    X = Dict{Int64, Vector{Float64}}(
+        1 => [0.0, 0.0],
+        2 => [1.0, 0.0],
+        3 => [1.0, 1.0],
+        4 => [0.0, 1.0])
+    update!([p1, p2], "geometry", 0.0 => X)
+    @test isapprox(p1("geometry", 0.0)[1], [0.0, 0.0])
+    @test isapprox(p2("geometry", 0.0)[1], [0.0, 0.0])
+    p1("geometry", 0.0)[1] = [1.0, 2.0]
+    @test isapprox(p2("geometry", 0.0)[1], [1.0, 2.0])
+end
+=#
+
+@testset "dict field depending from problems" begin
+    X = Dict{Int64, Vector{Float64}}(
+        1 => [0.0, 0.0],
+        2 => [1.0, 0.0],
+        3 => [1.0, 1.0],
+        4 => [0.0, 1.0],
+        5 => [0.0, 2.0],
+        6 => [1.0, 2.0],
+        7 => [1.0, 3.0],
+        8 => [0.0, 3.0])
+    p1 = Problem(Elasticity, "Body 1", 2)
+    p2 = Problem(Elasticity, "Body 2", 2)
+    e1 = Element(Quad4, 1, [1, 2, 3, 4])
+    e2 = Element(Quad4, 2, [5, 6, 7, 8])
+    push!(p1, e1)
+    push!(p2, e2)
+    update!(p1, "geometry", 0.0 => X)
+    update!(p2, "geometry", 0.0 => X)
+    @test isapprox(p1("geometry", 0.0)[1], [0.0, 0.0])
+    @test isapprox(p2("geometry", 0.0)[1], [0.0, 0.0])
+    p1("geometry", 0.0)[1] = [1.0, 2.0]
+    @test isapprox(p2("geometry", 0.0)[1], [1.0, 2.0])
+    @test isapprox(e1("geometry", 0.0)[1], [1.0, 2.0])
+end