diff --git a/src/JuliaFEM.jl b/src/JuliaFEM.jl
index 0697ce0..7bd7075 100644
--- a/src/JuliaFEM.jl
+++ b/src/JuliaFEM.jl
@@ -83,7 +83,7 @@ end
 include("assembly.jl")
 include("solver_utils.jl")
 include("solvers.jl")
-export AbstractSolver, Solver, Nonlinear, 
+export AbstractSolver, Solver, Nonlinear, NonlinearSolver, Linear, LinearSolver,
        get_unknown_field_name, get_formulation_type,
        get_field_problems, get_boundary_problems,
        get_field_assembly, get_boundary_assembly,
@@ -144,7 +144,7 @@ export get_mesh, get_model
 
 module Postprocess
 include("postprocess_utils.jl")
-export calc_nodal_values!, get_nodal_vector
+export calc_nodal_values!, get_nodal_vector, copy_field!
 include("postprocess_xdmf.jl")
 export XDMF, xdmf_new_result!, xdmf_save_field!, xdmf_save!
 end
diff --git a/src/assembly.jl b/src/assembly.jl
index 3497c19..6e8de19 100644
--- a/src/assembly.jl
+++ b/src/assembly.jl
@@ -14,15 +14,23 @@ type CAssembly
 end
 
 function optimize!(assembly::Assembly)
-    optimize!(assembly.mass_matrix)
-    optimize!(assembly.stiffness_matrix)
-    optimize!(assembly.force_vector)
+    optimize!(assembly.K)
+    optimize!(assembly.Kg)
+    optimize!(assembly.f)
+    optimize!(assembly.fg)
+    optimize!(assembly.C1)
+    optimize!(assembly.C2)
+    optimize!(assembly.D)
+    optimize!(assembly.g)
+    optimize!(assembly.c)
 end
 
 function append!(assembly::Assembly, sub_assembly::Assembly)
     append!(assembly.M, sub_assembly.M)
     append!(assembly.K, sub_assembly.K)
+    append!(assembly.Kg, sub_assembly.Kg)
     append!(assembly.f, sub_assembly.f)
+    append!(assembly.fg, sub_assembly.fg)
     append!(assembly.C1, sub_assembly.C1)
     append!(assembly.C2, sub_assembly.C2)
     append!(assembly.D, sub_assembly.D)
@@ -36,33 +44,38 @@ end
 function assemble_posthook!
 end
 
-function assemble!(problem::Problem, time::Real; empty_assembly::Bool=true)
+function assemble!(problem::Problem, time::Real)
+    if !isempty(problem.assembly)
+        warn("problem.assembly is not empty and assembling, are you sure you know what are you doing?")
+    end
     if method_exists(assemble_prehook!, Tuple{typeof(problem), Real})
         assemble_prehook!(problem, time)
     end
-    !problem.assembly.changed && return
-    empty_assembly && empty!(problem.assembly)
     for element in get_elements(problem)
         assemble!(problem.assembly, problem, element, time)
     end
-    problem.assembly.changed = true
     if method_exists(assemble_posthook!, Tuple{typeof(problem), Real})
         assemble_posthook!(problem, time)
     end
-    return
 end
 
-function assemble!(problem::Problem, time::Real, ::Type{Val{:mass_matrix}})
-    !isempty(problem.assembly.M) && return # assembly mass matrix only once
-    dim = get_unknown_field_dimension(problem)
+function assemble!(problem::Problem, time::Real, ::Type{Val{:mass_matrix}}; density=0.0, dual_basis=false, dim=0)
+    if !isempty(problem.assembly.M)
+        warn("problem.assembly.M is not empty and assembling, are you sure you know what are you doing?")
+    end
+    if dim == 0
+        dim = get_unknown_field_dimension(problem)
+    end
     for element in get_elements(problem)
-        haskey(element, "density") || error("Failed to assemble mass matrix, density not defined!")
+        if !haskey(element, "density") && density == 0.0
+            error("Failed to assemble mass matrix, density not defined!")
+        end
         nnodes = length(element)
         M = zeros(nnodes, nnodes)
         for ip in get_integration_points(element, 1)
             detJ = element(ip, time, Val{:detJ})
             N = element(ip, time)
-            rho = element("density", ip, time)
+            rho = haskey(element, "density") ? element("density", ip, time) : density
             M += ip.weight*rho*N'*N*detJ
         end
         gdofs = get_gdofs(problem, element)
diff --git a/src/contact.jl b/src/contact.jl
index 7a23288..10b0185 100644
--- a/src/contact.jl
+++ b/src/contact.jl
@@ -1,6 +1,15 @@
 # This file is a part of JuliaFEM.
 # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
 
+"""
+
+Parameters
+----------
+distval
+    a charasteristic measure to skip element pair, 0..5 => near, 10+ => far
+    5 means that distance of slave element midpoint and point to project
+    is 5 times larger than length of element
+"""
 type Contact <: BoundaryProblem
     dimension :: Int
     rotate_normals :: Bool
@@ -9,10 +18,15 @@ type Contact <: BoundaryProblem
     dual_basis :: Bool
     use_forwarddiff :: Bool
     minimum_active_set_size :: Int
+    distval :: Float64
+    store_fields :: Vector{ASCIIString}
 end
 
 function Contact()
-    return Contact(-1, false, false, false, true, false, 0)
+    default_fields = ["element area", "contact area", "weighted gap",
+        "contact pressure", "active nodes", "inactive nodes", "stick nodes",
+        "slip nodes", "complementarity condition", "contact error"]
+    return Contact(-1, false, false, false, true, false, 0, 5.0, default_fields)
 end
 
 function get_unknown_field_name(problem::Problem{Contact})
@@ -34,15 +48,14 @@ function assemble!(problem::Problem{Contact}, time::Real)
     dimension = Val{problem.properties.dimension}
     finite_sliding = Val{problem.properties.finite_sliding}
     friction = Val{problem.properties.friction}
-    dual_basis = Val{problem.properties.dual_basis}
     use_forwarddiff = Val{problem.properties.use_forwarddiff}
-    assemble!(problem, time, dimension, finite_sliding, friction, dual_basis, use_forwarddiff)
+    assemble!(problem, time, dimension, finite_sliding, friction, use_forwarddiff)
 end
 
-""" Frictionless 2d small sliding contact with dual basis without forwarddiff. """
-function assemble!(problem::Problem{Contact}, time::Real,
-                   ::Type{Val{1}}, ::Type{Val{false}}, ::Type{Val{false}},
-                   ::Type{Val{true}}, ::Type{Val{false}}; debug=false)
+""" Frictionless 2d small sliding contact without forwarddiff. """
+function assemble!(problem::Problem{Contact}, time::Float64,
+                   ::Type{Val{1}}, ::Type{Val{false}},
+                   ::Type{Val{false}}, ::Type{Val{false}}; debug=false)
 
     props = problem.properties
     field_dim = get_unknown_field_dimension(problem)
@@ -50,60 +63,80 @@ function assemble!(problem::Problem{Contact}, time::Real,
     slave_elements = get_slave_elements(problem)
 
     # 1. calculate nodal normals and tangents for slave element nodes j ∈ S
-    normals, tangents = calculate_normals(slave_elements, time, Val{1}; rotate_normals=props.rotate_normals)
-    update!(slave_elements, "normal", normals)
-    update!(slave_elements, "tangent", tangents)
+    normals, tangents = calculate_normals(slave_elements, time, Val{1};
+                        rotate_normals=props.rotate_normals)
+    update!(slave_elements, "normal", time => normals)
+    update!(slave_elements, "tangent", time => tangents)
 
     # 2. loop all slave elements
     for slave_element in slave_elements
 
+        nsl = length(slave_element)
         X1 = slave_element["geometry"](time)
         u1 = slave_element["displacement"](time)
         la1 = slave_element["reaction force"](time)
-        x1 = X1 + u1
         n1 = slave_element["normal"](time)
         t1 = slave_element["tangent"](time)
+        x1 = X1 + u1
         Q1_ = [n1[1] t1[1]]
         Q2_ = [n1[2] t1[2]]
         Z = zeros(2, 2)
         Q2 = [Q1_ Z; Z Q2_]
+        contact_area = 0.0
+        contact_error = 0.0
+
+        if "element area" in props.store_fields
+            element_area = 0.0
+            for ip in get_integration_points(slave_element, 3)
+                detJ = slave_element(ip, time, Val{:detJ})
+                w = ip.weight*detJ
+                element_area += w
+            end
+            update!(slave_element, "element area", time => element_area)
+        end
 
         # 3. loop all master elements
         for master_element in slave_element["master elements"](time)
 
+            nm = length(master_element)
             X2 = master_element["geometry"](time)
             u2 = master_element["displacement"](time)
             x2 = X2 + u2
 
+            norm(mean(X1) - X2[1]) / norm(X1[2] - X1[1]) < props.distval || continue
+            norm(mean(X1) - X2[2]) / norm(X1[2] - X1[1]) < props.distval || continue
+            
             # 3.1 calculate segmentation
             xi1a = project_from_master_to_slave(slave_element, X2[1], time)
-            xi1b = project_from_master_to_slave(slave_element, X2[end], time)
+            xi1b = project_from_master_to_slave(slave_element, X2[2], time)
             xi1 = clamp([xi1a; xi1b], -1.0, 1.0)
             l = 1/2*abs(xi1[2]-xi1[1])
             isapprox(l, 0.0) && continue # no contribution in this master element
 
             # 3.2. bi-orthogonal basis
-            nsl = length(slave_element)
-            nm = length(master_element)
             De = zeros(nsl, nsl)
             Me = zeros(nsl, nsl)
-            for ip in get_integration_points(slave_element, 3)
-                detJ = slave_element(ip, time, Val{:detJ})
-                w = ip.weight*detJ*l
-                xi = ip.coords[1]
-                xi_s = dot([1/2*(1-xi); 1/2*(1+xi)], xi1)
-                N1 = vec(get_basis(slave_element, xi_s, time))
-                De += w*diagm(N1)
-                Me += w*N1*N1'
+            Ae = zeros(nsl, nsl)
+            if props.dual_basis
+                for ip in get_integration_points(slave_element, 3)
+                    detJ = slave_element(ip, time, Val{:detJ})
+                    w = ip.weight*detJ*l
+                    xi = ip.coords[1]
+                    xi_s = dot([1/2*(1-xi); 1/2*(1+xi)], xi1)
+                    N1 = vec(get_basis(slave_element, xi_s, time))
+                    De += w*diagm(N1)
+                    Me += w*N1*N1'
+                end
+                Ae = De*inv(Me)
+            else
+                Ae = eye(nsl)
             end
-            Ae = De*inv(Me)
             
             # 3.3. loop integration points of one integration segment and calculate
             # local mortar matrices
             fill!(De, 0.0)
             fill!(Me, 0.0)
             ge = zeros(field_dim*nsl)
-            lae = zeros(field_dim*nsl)
             for ip in get_integration_points(slave_element, 3)
                 detJ = slave_element(ip, time, Val{:detJ})
                 w = ip.weight*detJ*l
@@ -119,11 +152,14 @@ function assemble!(problem::Problem{Contact}, time::Real,
                 X_m = N2*X2 
                 De += w*Phi*N1'
                 Me += w*Phi*N2'
-                x_s = X_s + N1*u1
-                x_m = X_m + N2*u2
+                u_s = N1*u1
+                u_m = N2*u2
+                x_s = X_s + u_s
+                x_m = X_m + u_m
                 la_s = Phi*la1
                 ge += w*vec((x_m-x_s)*Phi')
-                lae += w*vec(la_s*Phi')
+                contact_area += w
+                contact_error += 1/2*w*dot(n_s, x_s-x_m)^2
             end
 
             # add contribution to contact virtual work
@@ -139,63 +175,117 @@ function assemble!(problem::Problem{Contact}, time::Real,
             end
             add!(problem.assembly.C1, sdofs, sdofs, D2)
             add!(problem.assembly.C1, sdofs, mdofs, -M2)
-
             add!(problem.assembly.C2, sdofs, sdofs, Q2'*D2)
             add!(problem.assembly.C2, sdofs, mdofs, -Q2'*M2)
             add!(problem.assembly.g, sdofs, Q2'*ge)
-            add!(problem.assembly.c, sdofs, Q2'*lae)
 
         end # master elements done
 
+        if "contact area" in props.store_fields
+            update!(slave_element, "contact area", time => contact_area)
+        end
+
+        if "contact error" in props.store_fields
+            update!(slave_element, "contact error", time => contact_error)
+        end
+
     end # slave elements done, contact virtual work ready
 
     S = sort(collect(keys(normals))) # slave element nodes
+    weighted_gap = Dict{Int64, Vector{Float64}}()
+    contact_pressure = Dict{Int64, Vector{Float64}}()
+    complementarity_condition = Dict{Int64, Vector{Float64}}()
+    is_active = Dict{Int64, Int}()
+    is_inactive = Dict{Int64, Int}()
+    is_slip = Dict{Int64, Int}()
+    is_stick = Dict{Int64, Int}()
+
+    g = full(problem.assembly.g)
+    la = problem.assembly.la
+
+    # active / inactive node detection
+    for j in S
+        dofs = [2*(j-1)+1, 2*(j-1)+2]
+        weighted_gap[j] = g[dofs]
+        if length(la) != 0
+            p = dot(normals[j], la[dofs])
+            t = dot(tangents[j], la[dofs])
+            contact_pressure[j] = [p, t]
+        else
+            contact_pressure[j] = [0.0, 0.0]
+        end
+        complementarity_condition[j] = contact_pressure[j] - weighted_gap[j]
+        if complementarity_condition[j][1] < 0
+            is_inactive[j] = 1
+            is_active[j] = 0
+            is_slip[j] = 0
+            is_stick[j] = 0
+        else
+            is_inactive[j] = 0
+            is_active[j] = 1
+            is_slip[j] = 1
+            is_stick[j] = 0
+        end
+    end
+
+    if "weighted gap" in props.store_fields
+        update!(slave_elements, "weighted gap", time => weighted_gap)
+    end
+    if "contact pressure" in props.store_fields
+        update!(slave_elements, "contact pressure", time => contact_pressure)
+    end
+    if "complementarity condition" in props.store_fields
+        update!(slave_elements, "complementarity condition", time => complementarity_condition)
+    end
+    if "active nodes" in props.store_fields
+        update!(slave_elements, "active nodes", time => is_active)
+    end
+    if "inactive nodes" in props.store_fields
+        update!(slave_elements, "inactive nodes", time => is_inactive)
+    end
+    if "stick nodes" in props.store_fields
+        update!(slave_elements, "stick nodes", time => is_stick)
+    end
+    if "slip nodes" in props.store_fields
+        update!(slave_elements, "slip nodes", time => is_slip)
+    end
+
+    info("# | active | inactive | stick | slip | gap | pres | comp")
+    for j in S
+        str1 = "$j | $(is_active[j]) | $(is_inactive[j]) | $(is_stick[j]) | $(is_slip[j]) | "
+        str2 = "$(round(weighted_gap[j], 3)) | $(round(contact_pressure[j], 3)) | $(round(complementarity_condition[j], 3))"
+        info(str1 * str2)
+    end
+
+    # solve variational inequality
+    
     C1 = sparse(problem.assembly.C1)
     ndofs = size(C1, 1)
-    debug && info("ndofs = $ndofs")
     C2 = sparse(problem.assembly.C2)
     D = spzeros(ndofs, ndofs)
-    g = sparse(problem.assembly.g)
-    g = full(g)
-    c = sparse(problem.assembly.c)
-    c = full(c)
-    debug && info("Contact slave nodes: $S")
 
     # constitutive modelling in tangent direction, frictionless contact
     for j in S
         dofs = [2*(j-1)+1, 2*(j-1)+2]
-        C2[dofs[2],:] = 0.0
-        g[dofs[2]] = 0.0
-        D[dofs[2], dofs] = tangents[j]
+        if (is_active[j] == 1) && (is_slip[j] == 1)
+            info("$j is in active/slip, removing tangential constraint $(dofs[2])")
+            C2[dofs[2],:] = 0.0
+            g[dofs[2]] = 0.0
+            D[dofs[2], dofs] = tangents[j]
+        end
     end
-    debug && info("Constitutive modelling ready")
 
-    # active / inactive node detection
-    A = Set()
-    I = Set()
-    la = problem.assembly.la
+    # remove inactive nodes from assembly
     for j in S
         dofs = [2*(j-1)+1, 2*(j-1)+2]
-        
-        Cn = -g[dofs[1]]
-        if length(la) != 0
-            Cn += dot(normals[j], la[dofs])
-            debug && info("slave $j: $(normals[j]) | $(la[dofs]) | $(c[dofs]) | $(g[dofs]) | $Cn")
-        else
-            debug && info("slave $j: $(normals[j]) |             | $(c[dofs]) | $(g[dofs]) | $Cn")
-        end
-        if Cn < 0
-            push!(I, j)
-            debug && info("slave $j INACTIVE")
+        if is_inactive[j] == 1
+            info("$j is inactive, removing dofs $dofs")
             C1[dofs,:] = 0.0
             C2[dofs,:] = 0.0
             D[dofs,:] = 0.0
             g[dofs,:] = 0.0
-        else
-            push!(A, j)
         end
     end
-    debug && info("active nodes: $A, inactive nodes: $I")
 
     problem.assembly.C1 = C1
     problem.assembly.C2 = C2
diff --git a/src/contact_2d_autodiff.jl b/src/contact_2d_autodiff.jl
new file mode 100644
index 0000000..7849b73
--- /dev/null
+++ b/src/contact_2d_autodiff.jl
@@ -0,0 +1,272 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
+
+""" Find segment from slave element corresponding to master element nodes.
+
+Parameters
+----------
+x1_, n1_
+    slave element geometry and normal direction
+x2
+    master element node to project onto slave
+
+Returns
+-------
+xi
+    dimensionless coordinate on slave corresponding to
+    projected master
+ 
+"""
+function project_from_master_to_slave{E<:MortarElements2D}(
+    slave_element::Element{E}, x1_::DVTI, n1_::DVTI, x2::Vector;
+    tol=1.0e-10, max_iterations=20)
+
+    x1(xi1) = vec(get_basis(E, xi1))*x1_
+    dx1(xi1) = vec(get_dbasis(E, xi1))*x1_
+    n1(xi1) = vec(get_basis(E, xi1))*n1_
+    dn1(xi1) = vec(get_dbasis(E, xi1))*n1_
+    cross2(a, b) = cross([a; 0], [b; 0])[3]
+    R(xi1) = cross2(x1(xi1)-x2, n1(xi1))
+    dR(xi1) = cross2(dx1(xi1), n1(xi1)) + cross2(x1(xi1)-x2, dn1(xi1))
+
+    xi1 = 0.0
+    dxi1 = 0.0
+    for i=1:max_iterations
+        dxi1 = -R(xi1)/dR(xi1)
+        xi1 += dxi1
+        if norm(dxi1) < tol
+            return xi1
+        end
+    end
+
+    info("x1 = $(ForwardDiff.get_value(x1_.data))")
+    info("n1 = $(ForwardDiff.get_value(n1_.data))")
+    info("x2 = $(ForwardDiff.get_value(x2))")
+    info("xi1 = $(ForwardDiff.get_value(xi1)), dxi1 = $(ForwardDiff.get_value(dxi1))")
+    info("-R(xi1) = $(ForwardDiff.get_value(-R(xi1)))")
+    info("dR(xi1) = $(ForwardDiff.get_value(dR(xi1)))")
+    error("find projection from master to slave: did not converge")
+
+end
+
+function project_from_slave_to_master{E<:MortarElements2D}(
+    master_element::Element{E}, x1::Vector, n1::Vector, x2_::DVTI;
+    tol=1.0e-10, max_iterations=20)
+
+    x2(xi2) = vec(get_basis(E, xi2))*x2_
+    dx2(xi2) = vec(get_dbasis(E, xi2))*x2_
+    cross2(a, b) = cross([a; 0], [b; 0])[3]
+    R(xi2) = cross2(x2(xi2)-x1, n1)
+    dR(xi2) = cross2(dx2(xi2), n1)
+
+    xi2 = 0.0
+    dxi2 = 0.0
+    for i=1:max_iterations
+        dxi2 = -R(xi2) / dR(xi2)
+        xi2 += dxi2
+        if norm(dxi2) < tol
+            return xi2
+        end
+    end
+
+    error("find projection from slave to master: did not converge, last val: $xi2 and $dxi2")
+
+end
+
+""" Assemble Mortar problem for two-dimensional problems, i.e. for Seg2 and Seg3 elements. """
+function assemble!(problem::Problem{Mortar}, time::Real, ::Type{Val{2}})
+
+    props = problem.properties
+    field_dim = get_unknown_field_dimension(problem)
+    field_name = get_parent_field_name(problem)
+
+    function calculate_interface(x::Vector)
+
+        ndofs = round(Int, length(x)/2)
+        nnodes = round(Int, ndofs/field_dim)
+        u = reshape(x[1:ndofs], field_dim, nnodes)
+        la = reshape(x[ndofs+1:end], field_dim, nnodes)
+        fc = zeros(u)
+        gap = zeros(u)
+        C = zeros(la)
+        S = Set{Int64}()
+
+        # 1. update nodal normals for slave elements
+        Q = [0.0 -1.0; 1.0 0.0]
+        normals = zeros(u)
+        for element in get_elements(problem)
+            haskey(element, "master elements") || continue
+            conn = get_connectivity(element)
+            push!(S, conn...)
+            gdofs = get_gdofs(element, field_dim)
+            X_el = element("geometry", time)
+            u_el = Field(Vector[u[:,i] for i in conn])
+            x_el = X_el + u_el
+            for ip in get_integration_points(element, Val{3})
+                dN = get_dbasis(element, ip)
+                N = element(ip, time)
+                t = sum([kron(dN[:,i], x_el[i]') for i=1:length(x_el)])
+                normals[:, conn] += ip.weight*Q*t'*N
+            end
+        end
+        for i in 1:size(normals,2)
+            normals[:,i] /= norm(normals[:,i])
+        end
+        # swap element normals in 2d if they point to inside of body
+        if props.rotate_normals
+            for i=1:size(normals,2)
+                normals[:,i] = -normals[:,i]
+            end
+        end
+
+        # 2. loop all slave elements
+        for slave_element in get_elements(problem)
+            haskey(slave_element, "master elements") || continue
+
+            slave_element_nodes = get_connectivity(slave_element)
+            X1 = slave_element("geometry", time)
+            u1 = Field(Vector[u[:,i] for i in slave_element_nodes])
+            x1 = X1 + u1
+            la1 = Field(Vector[la[:,i] for i in slave_element_nodes])
+            n1 = Field(Vector[normals[:,i] for i in slave_element_nodes])
+            nnodes = size(slave_element, 2)
+            update!(slave_element, "normals", time => ForwardDiff.get_value(n1.data))
+
+            # 3. loop all master elements
+            for master_element in slave_element["master elements"]
+
+                master_element_nodes = get_connectivity(master_element)
+                X2 = master_element("geometry", time)
+                u2 = Field(Vector[u[:,i] for i in master_element_nodes])
+                x2 = X2 + u2
+
+                x1_midpoint = 1/2*(x1[1]+x1[2])
+                x2_midpoint = 1/2*(x2[1]+x2[2])
+                distance = ForwardDiff.get_value(norm(x2_midpoint - x1_midpoint))
+                distance > props.maximum_distance && continue
+
+                # calculate segmentation: we care only about endpoints
+                # note: these are quadratic/cubic functions, analytical solution possible
+                xi1a = -Inf
+                xi1b = -Inf
+                try
+                    xi1a = project_from_master_to_slave(slave_element, x1, n1, x2[1])
+                    xi1b = project_from_master_to_slave(slave_element, x1, n1, x2[end])
+                catch
+                    info("failed to create projection!!!!")
+                    # TODO
+                    continue
+                end
+                xi1 = clamp([xi1a; xi1b], -1.0, 1.0)
+                l = 1/2*abs(xi1[2]-xi1[1])
+                isapprox(l, 0.0) && continue # no contribution in this master element
+
+                De = zeros(nnodes, nnodes)
+                Me = zeros(nnodes, nnodes)
+                for ip in get_integration_points(slave_element, Val{5})
+                    # jacobian of slave element in deformed state
+                    dN = get_dbasis(slave_element, ip)
+                    j = sum([kron(dN[:,i], x1[i]') for i=1:length(x1)])
+                    w = ip.weight*norm(j)*l
+                    xi_s = dot([1/2*(1-ip.xi); 1/2*(1+ip.xi)], xi1)
+                    N1 = get_basis(slave_element, xi_s)
+                    De += w*diagm(vec(N1))
+                    Me += w*N1'*N1
+                end
+                Ae = De*inv(Me)
+
+                slave_dofs = get_gdofs(slave_element, field_dim)
+                master_dofs = get_gdofs(master_element, field_dim)
+
+                # 4. loop integration points of segment
+                for ip in get_integration_points(slave_element, Val{5})
+                    # jacobian of slave element in deformed state
+                    dN = get_dbasis(slave_element, ip)
+                    j = sum([kron(dN[:,i], x1[i]') for i=1:length(x1)])
+                    w = ip.weight*norm(j)*l
+
+                    # project gauss point from slave element to master element
+                    xi_s = dot([1/2*(1-ip.xi); 1/2*(1+ip.xi)], xi1)
+                    N1 = vec(get_basis(slave_element, xi_s))
+                    x_s = N1*x1 # coordinate in gauss point
+                    n_s = N1*n1 # normal direction in gauss point
+                    t_s = Q'*n_s # tangent direction in gauss point
+                    xi_m = project_from_slave_to_master(master_element, x_s, n_s, x2)
+                    N2 = vec(get_basis(master_element, xi_m))
+                    x_m = N2*x2
+                    Phi = Ae*N1
+
+                    la_s = Phi*la1 # traction force in gauss point
+                    gn = props.gap_sign*dot(n_s, x_s - x_m) # normal gap
+
+                    fc[:,slave_element_nodes] += w*la_s*N1'
+                    fc[:,master_element_nodes] -= w*la_s*N2'
+                    gap[1,slave_element_nodes] += w*gn*Phi'
+                    #gap[1,slave_element_nodes] += w*gn*N1'
+
+                end # done integrating segment
+
+            end # master elements done
+
+        end # slave elements done
+
+        # at this point we have calculated contact force fc and gap for all slave elements.
+        # next task is to find out are they in contact or not and remove inactive nodes
+
+        nzgap = sort(nonzeros(sparse(ForwardDiff.get_value(gap))))
+        info("gap: $nzgap")
+
+        for (i, j) in enumerate(sort(collect(S)))
+            if j in props.always_inactive
+                info("special node $j always inactive")
+                C[:,j] = la[:,j]
+                continue
+            end
+            n = normals[:,j]
+            t = Q'*n
+            lan = dot(n, la[:,j])
+            lat = dot(t, la[:,j])
+
+            if lan - gap[1, j] > 0
+                info("set node $j active, normal direction = $(ForwardDiff.get_value(n)), tangent plane = $(ForwardDiff.get_value(t))")
+                C[1,j] += gap[1, j]
+                C[2,j] += lat
+            else
+                C[:,j] = la[:,j]
+            end
+        end
+
+        return vec([fc C])
+
+    end
+
+    # x doesn't mean deformed configuration here
+    x = [problem.assembly.u; problem.assembly.la]
+    ndofs = round(Int, length(x)/2)
+    A, allresults = ForwardDiff.jacobian(calculate_interface, x,
+                            ForwardDiff.AllResults, cache=autodiffcache)
+    b = -ForwardDiff.value(allresults)
+
+    A = sparse(A)
+    b = sparse(b)
+    SparseMatrix.droptol!(A, 1.0e-12)
+    SparseMatrix.droptol!(b, 1.0e-12)
+
+    K = A[1:ndofs,1:ndofs]
+    C1 = transpose(A[1:ndofs,ndofs+1:end])
+    C2 = A[ndofs+1:end,1:ndofs]
+    D = A[ndofs+1:end,ndofs+1:end]
+    f = b[1:ndofs]
+    g = b[ndofs+1:end]
+
+    empty!(problem.assembly)
+    add!(problem.assembly.K, K)
+    add!(problem.assembly.C1, C1)
+    add!(problem.assembly.C2, C2)
+    add!(problem.assembly.D, D)
+    add!(problem.assembly.f, f)
+    add!(problem.assembly.g, g)
+
+    return problem.assembly
+
+end
diff --git a/src/elasticity.jl b/src/elasticity.jl
index 9fbd435..81f10ba 100644
--- a/src/elasticity.jl
+++ b/src/elasticity.jl
@@ -10,10 +10,11 @@ type Elasticity <: FieldProblem
     formulation :: Symbol
     finite_strain :: Bool
     geometric_stiffness :: Bool
+    store_fields :: Vector{ASCIIString}
 end
 function Elasticity()
     # formulations: plane_stress, plane_strain, continuum
-    return Elasticity(:continuum, false, false)
+    return Elasticity(:continuum, false, false, [])
 end
 
 function get_unknown_field_name(problem::Problem{Elasticity})
@@ -84,7 +85,6 @@ function assemble{El<:Union{Tri3,Tri6,Quad4}}(problem::Problem{Elasticity}, elem
         end
 
         strain_vec = [strain[1,1]; strain[2,2]; strain[1,2]]
-        update!(ip, "strain", time => strain_vec)
 
         # calculate stress
         E = element("youngs modulus", ip, time)
@@ -104,7 +104,12 @@ function assemble{El<:Union{Tri3,Tri6,Quad4}}(problem::Problem{Elasticity}, elem
         end
         # calculate stress
         stress_vec = D * ([1.0, 1.0, 2.0] .* strain_vec)
-        update!(ip, "stress", time => stress_vec)
+
+        "strain" in props.store_fields && update!(ip, "strain", time => strain_vec)
+        "stress" in props.store_fields && update!(ip, "stress", time => stress_vec)
+        "stress 11" in props.store_fields && update!(ip, "stress 11", time => stress_vec[1])
+        "stress 22" in props.store_fields && update!(ip, "stress 22", time => stress_vec[2])
+        "stress 12" in props.store_fields && update!(ip, "stress 12", time => stress_vec[3])
 
         Km += w*BL'*D*BL
 
@@ -400,7 +405,6 @@ function assemble{El<:Union{Tet4, Tet10, Hex8}}(problem::Problem{Elasticity}, el
         end
 
         strain_vec = [strain[1,1]; strain[2,2]; strain[3,3]; strain[1,2]; strain[2,3]; strain[1,3]]
-        update!(ip, "strain", time => strain_vec)
 
         # calculate stress
         E = element("youngs modulus", ip, time)
@@ -413,7 +417,15 @@ function assemble{El<:Union{Tet4, Tet10, Hex8}}(problem::Problem{Elasticity}, el
             0.0 0.0 0.0 0.0 0.5-nu 0.0
             0.0 0.0 0.0 0.0 0.0 0.5-nu]
         stress_vec = D * ([1.0, 1.0, 1.0, 2.0, 2.0, 2.0].*strain_vec)
-        update!(ip, "stress", time => stress_vec)
+
+        "strain" in props.store_fields && update!(ip, "strain", time => strain_vec)
+        "stress" in props.store_fields && update!(ip, "stress", time => stress_vec)
+        "stress 11" in props.store_fields && update!(ip, "stress 11", time => stress_vec[1])
+        "stress 22" in props.store_fields && update!(ip, "stress 22", time => stress_vec[2])
+        "stress 33" in props.store_fields && update!(ip, "stress 33", time => stress_vec[3])
+        "stress 12" in props.store_fields && update!(ip, "stress 12", time => stress_vec[4])
+        "stress 23" in props.store_fields && update!(ip, "stress 23", time => stress_vec[5])
+        "stress 13" in props.store_fields && update!(ip, "stress 13", time => stress_vec[6])
 
         Km += w*BL'*D*BL
 
diff --git a/src/elements.jl b/src/elements.jl
index 49a04d8..8191ee5 100644
--- a/src/elements.jl
+++ b/src/elements.jl
@@ -18,7 +18,11 @@ function Element{E<:AbstractElement}(::Type{E}, connectivity=[], integration_poi
 end
 
 function getindex(element::Element, field_name::ASCIIString)
-    element.fields[field_name]
+    return element.fields[field_name]
+end
+
+function setindex!(element::Element, data::Field, field_name::ASCIIString)
+    element.fields[field_name] = data
 end
 
 function setindex!(element::Element, data, field_name::ASCIIString)
@@ -30,7 +34,7 @@ function call(element::Element, field_name::ASCIIString, time)
 end
 
 function call(element::Element, ip, time)
-    get_basis(element, ip, time)
+    return get_basis(element, ip, time)
 end
 
 function call(element::Element, ip, time, ::Type{Val{:Jacobian}})
@@ -60,44 +64,37 @@ function call(element::Element, ip, time, ::Type{Val{:Grad}})
 end
 
 function call(element::Element, field_name::ASCIIString, ip, time, ::Type{Val{:Grad}})
-    element(ip, time, Val{:Grad})*element[field_name](time)
+    return element(ip, time, Val{:Grad})*element[field_name](time)
 end
 
 function call(element::Element, field_name::ASCIIString, time)
     return element[field_name](time)
 end
 
-function call(element::Element, field_name::ASCIIString, ip, time::Real)
-    field = element(field_name, time)
-    isa(field, DCTI) && return field.data
-    basis = element(ip, time)
-    n = length(element)
-    m = length(field)
-    @assert n == m
-    return sum([field[i]*basis[i] for i=1:n])
+function call(element::Element, field_name::ASCIIString, ip, time::Float64)
+    field = element[field_name]
+    return call(element, field, ip, time)
 end
 
-#function get_jacobian{E}(element::Element{E}, xi::Vector, time=0.0)
-#    element(xi, time, Val{:Jacobian})
-#end
-#function get_basis(element::Element, xi::Vector, time=0.0)
-#    get_basis(element.properties, xi, time)
-#end
-#function get_dbasis(element::Element, xi::Vector, time=0.0)
-#    get_dbasis(element.properties, xi, time)
-#end
-#function get_integration_points{E}(element::Element{E})
-#    get_integration_points(element.properties)
-#end
-#function length{E}(element::Element{E})
-#    length(element.properties)
-#end
-#function size{E}(element::Element{E})
-#    size(element.properties)
-#end
+function call(element::Element, field::DCTI, ip, time::Float64)
+    return field.data
+end
+
+function call(element::Element, field::CVTV, ip, time::Float64)
+    return field(ip, time)
+end
+
+function call(element::Element, field::Field, ip, time::Float64)
+    field_ = field(time)
+    basis = element(ip, time)
+    n = length(element)
+    m = length(field_)
+    @assert n == m
+    return sum([field_[i]*basis[i] for i=1:n])
+end
 
 function size(element::Element, dim::Int)
-    size(element)[dim]
+    return size(element)[dim]
 end
 
 """ Update element field based on a dictionary of nodal data and connectivity information.
@@ -115,10 +112,64 @@ function update!(element::Element, field_name::ASCIIString, data::Dict)
     element[field_name] = [data[i] for i in get_connectivity(element)]
 end
 
-function update!(element::Element, field_name::ASCIIString, datas::Union{Real, Vector, Pair}...)
+function update!{K,V}(element::Element, field_name::ASCIIString, 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)
+end
+
+function update!(element::Element, field_name::ASCIIString, datas::Union{Real, Vector, Pair{Float64, Union{Real, Vector{Any}}}}...)
     for data in datas
         if haskey(element, field_name)
             update!(element[field_name], data)
+        else
+            if length(data) != length(element)
+                update!(element, field_name, DCTI(data))
+            else
+                element[field_name] = data
+            end
+        end
+    end
+end
+
+function update!(element::Element, field_name::ASCIIString, data::Pair{Float64, Vector{Any}})
+    if haskey(element, field_name)
+        update!(element[field_name], data)
+    else
+        element[field_name] = data
+    end
+end
+
+function update!(element::Element, field_name::ASCIIString, data::Pair{Float64, Vector{Int64}})
+    if haskey(element, field_name)
+        update!(element[field_name], data)
+    else
+        element[field_name] = data
+    end
+end
+
+function update!(element::Element, field_name::ASCIIString, data::Pair{Float64, Vector{Vector{Float64}}})
+    if haskey(element, field_name)
+        update!(element[field_name], data)
+    else
+        element[field_name] = data
+    end
+end
+
+function update!(element::Element, field_name::ASCIIString, data::Pair{Float64, Float64})
+    if haskey(element, field_name)
+        update!(element[field_name], data)
+    else
+        element[field_name] = data
+    end
+end
+
+function update!(element::Element, field_name::ASCIIString, data::Union{Float64, Vector})
+    if haskey(element, field_name)
+        update!(element[field_name], data)
+    else
+        if length(data) != length(element)
+            update!(element, field_name, DCTI(data))
         else
             element[field_name] = data
         end
@@ -135,6 +186,14 @@ function update!(element::Element, datas::Pair...)
     end
 end
 
+function update!(element::Element, field_name::ASCIIString, data::Function)
+    element[field_name] = data
+end
+
+function update!(element::Element, field_name::ASCIIString, field::Field)
+    element[field_name] = field
+end
+
 function update!(elements::Vector, field_name::ASCIIString, data)
     for element in elements
         update!(element, field_name, data)
diff --git a/src/fields.jl b/src/fields.jl
index 38fb666..6839f26 100644
--- a/src/fields.jl
+++ b/src/fields.jl
@@ -1,8 +1,6 @@
 # This file is a part of JuliaFEM.
 # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
 
-# https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/notebooks/2015-06-14-data-structures.ipynb
-
 abstract AbstractField
 
 abstract Discrete <: AbstractField
@@ -12,7 +10,6 @@ abstract Variable <: AbstractField
 abstract TimeVariant <: AbstractField
 abstract TimeInvariant <: AbstractField
 
-
 type Field{A<:Union{Discrete,Continuous}, B<:Union{Constant,Variable}, C<:Union{TimeVariant,TimeInvariant}}
     data
 end
@@ -49,11 +46,11 @@ type Basis
     dbasis :: Function
 end
 
-function Base.call(basis::Basis, xi::Vector)
+function call(basis::Basis, xi::Vector)
     basis.basis(xi)
 end
 
-function Base.call(basis::Basis, xi::Vector, ::Type{Val{:grad}})
+function call(basis::Basis, xi::Vector, ::Type{Val{:grad}})
     basis.dbasis(xi)
 end
 
@@ -102,7 +99,7 @@ function Field{T}(data::Pair{Float64, Vector{T}}...)
     return DVTV([Increment{Vector{T}}(d[1], d[2]) for d in data])
 end
 
-function Base.convert{T}(::Type{DCTV}, data::Pair{Real, Vector{T}}...)
+function convert{T}(::Type{DCTV}, data::Pair{Real, Vector{T}}...)
     return DCTV([Increment{Vector{T}}(d[1], d[2]) for d in data])
 end
 
@@ -114,7 +111,7 @@ julia> t0 = 0.0; t1=1.0; y0 = 0.0; y1 = 1.0
 julia> f = DCTV(t0 => y0, t1 => y1)
 
 """
-function Base.convert{T,v<:Real}(::Type{DCTV}, data::Pair{v, T}...)
+function convert{T,v<:Real}(::Type{DCTV}, data::Pair{v, T}...)
     return DCTV([Increment(d[1],d[2]) for d in data])
 end
 #function Base.convert(::Type{DCTV}, data::Pair{Real, Any}...)
@@ -234,16 +231,16 @@ function Base.(:*)(T::Vector, f::DVTI)
     return sum([T[i]*f[i] for i=1:length(f)])
 end
 
-function Base.vec(field::DVTI)
+function vec(field::DVTI)
     return [field.data...;]
 end
 
-function Base.vec(field::DCTV)
+function vec(field::DCTV)
     info("trying to vectorize $field")
     error("does not make sense")
 end
 
-function Base.endof(field::Field)
+function endof(field::Field)
     return endof(field.data)
 end
 
@@ -251,22 +248,22 @@ end
 #    return Increment(reshape(data, round(Int, length(data)/length(increment)), length(increment)))
 #end
 
-function Base.similar{T}(field::DVTI, data::Vector{T})
+function similar{T}(field::DVTI, data::Vector{T})
     n = length(field.data)
     data = reshape(data, round(Int, length(data)/n), n)
     newdata = Vector[data[:,i] for i=1:n]
     return typeof(field)(newdata)
 end
 
-function Base.start(::DVTI)
+function start(::DVTI)
     return 1
 end
 
-function Base.next(f::DVTI, state)
+function next(f::DVTI, state)
     return f.data[state], state+1
 end
 
-function Base.done(f::DVTI, s)
+function done(f::DVTI, s)
     return s > length(f.data)
 end
 
@@ -301,22 +298,26 @@ end
 
 ### Accessing continuous fields
 
-function Base.call(field::CVTI, xi::Vector)
-    field.data(xi)
+function call(field::CVTI, xi::Vector)
+    return field.data(xi)
 end
 
-function Base.call(field::CVTI, xi::Vector, ::Type{Val{:grad}})
-    field.data(xi, Val{:grad})
+function call(field::CVTV, xi, time::Float64)
+    return field.data(xi, time)
 end
 
-function Base.convert(::Type{Basis}, field::CVTI)
-    return field.data
+function call(field::CVTI, xi::Vector, ::Type{Val{:Grad}})
+    return field.data(xi, Val{:Grad})
 end
 
-function Base.call(field::CCTV, time::Number)
+function call(field::CCTV, time::Float64)
     return field.data(time)
 end
 
+function convert(::Type{Basis}, field::CVTI)
+    return field.data
+end
+
 ### Interpolation
 
 """ Interpolate time-invariant field in time direction. """
diff --git a/src/mortar.jl b/src/mortar.jl
index 67b9cc8..b6787ea 100644
--- a/src/mortar.jl
+++ b/src/mortar.jl
@@ -5,12 +5,15 @@ type Mortar <: BoundaryProblem
     dimension :: Int
     rotate_normals :: Bool
     adjust :: Bool
-    tolerance :: Float64
     dual_basis :: Bool
+    use_forwarddiff :: Bool
+    distval :: Float64
+    store_fields :: Vector{ASCIIString}
 end
 
 function Mortar()
-    return Mortar(-1, false, false, 0.0, false)
+    default_fields = []
+    return Mortar(-1, false, false, false, false, Inf, default_fields)
 end
 
 function get_unknown_field_name(problem::Problem{Mortar})
@@ -19,6 +22,13 @@ end
 
 function get_formulation_type(problem::Problem{Mortar})
     return :incremental
+    #=
+    if problem.properties.use_forwarddiff
+        return :forwarddiff
+    else
+        return :incremental
+    end
+    =#
 end
 
 typealias MortarElements2D Union{Seg2, Seg3}
@@ -52,7 +62,25 @@ function project_from_master_to_slave{E<:MortarElements2D}(slave_element::Elemen
     dn1(xi1) = vec(get_dbasis(slave_element, [xi1], time))*n1_
     R(xi1) = cross2(x1(xi1)-x2, n1(xi1))
     dR(xi1) = cross2(dx1(xi1), n1(xi1)) + cross2(x1(xi1)-x2, dn1(xi1))
-    xi1 = newton(R, dR, 0.0)
+    xi1 = nothing
+    try
+        xi1 = newton(R, dR, 0.0)
+    catch
+        warn("projection from master to slave failed with following arguments:")
+        warn("slave element x1: $x1_")
+        warn("slave element n1: $n1_")
+        warn("master element x2: $x2")
+        warn("time: $time")
+        len = norm(x1_[2] - x1_[1])
+        midpnt = mean(x1_)
+        dist = norm(midpnt - x2)
+        distval = dist/len
+        warn("midpoint of slave element: $midpnt")
+        warn("length of slave element: $len")
+        warn("distance between midpoint of slave element and x2: $dist")
+        warn("charasteristic measure: $distval")
+        rethrow()
+    end
     return xi1
 end
 
@@ -109,16 +137,18 @@ function calculate_normals!(elements, time, ::Type{Val{1}}; rotate_normals=false
     end
 end
 
-function assemble!(problem::Problem{Mortar}, time::Real)
+function assemble!(problem::Problem{Mortar}, time::Float64)
     if problem.properties.dimension == -1
         problem.properties.dimension = dim = size(first(problem.elements), 1)
         info("assuming dimension of mesh tie surface is $dim")
         info("if this is wrong set is manually using problem.properties.dimension")
     end
-    assemble!(problem, time, Val{problem.properties.dimension})
+    dimension = Val{problem.properties.dimension}
+    use_forwarddiff = Val{problem.properties.use_forwarddiff}
+    assemble!(problem, time, dimension, use_forwarddiff)
 end
 
-function assemble!(problem::Problem{Mortar}, time::Real, ::Type{Val{1}})
+function assemble!(problem::Problem{Mortar}, time::Float64, ::Type{Val{1}}, ::Type{Val{false}})
 
     props = problem.properties
     field_dim = get_unknown_field_dimension(problem)
@@ -129,32 +159,51 @@ function assemble!(problem::Problem{Mortar}, time::Real, ::Type{Val{1}})
     normals, tangents = calculate_normals(slave_elements, time, Val{1};
                                           rotate_normals=props.rotate_normals)
     update!(slave_elements, "normal", normals)
+    update!(slave_elements, "tangent", tangents)
 
     # 2. loop all slave elements
     for slave_element in slave_elements
-        haskey(slave_element, "master elements") || continue
 
+        nsl = length(slave_element)
         X1 = slave_element["geometry"](time)
         n1 = slave_element["normal"](time)
 
         # 3. loop all master elements
         for master_element in slave_element["master elements"](time)
 
+            nm = length(master_element)
             X2 = master_element["geometry"](time)
 
             # 3.1 calculate segmentation
             xi1a = project_from_master_to_slave(slave_element, X2[1], time)
-            xi1b = project_from_master_to_slave(slave_element, X2[end], time)
+            xi1b = project_from_master_to_slave(slave_element, X2[2], time)
             xi1 = clamp([xi1a; xi1b], -1.0, 1.0)
             l = 1/2*abs(xi1[2]-xi1[1])
             isapprox(l, 0.0) && continue # no contribution in this master element
 
+            # 3.2. bi-orthogonal basis
+            De = zeros(nsl, nsl)
+            Me = zeros(nsl, nsl)
+            Ae = zeros(nsl, nsl)
+            if props.dual_basis
+                for ip in get_integration_points(slave_element, 3)
+                    detJ = slave_element(ip, time, Val{:detJ})
+                    w = ip.weight*detJ*l
+                    xi = ip.coords[1]
+                    xi_s = dot([1/2*(1-xi); 1/2*(1+xi)], xi1)
+                    N1 = vec(get_basis(slave_element, xi_s, time))
+                    De += w*diagm(N1)
+                    Me += w*N1*N1'
+                end
+                Ae = De*inv(Me)
+            else
+                Ae = eye(nsl)
+            end
+
             # 3.3. loop integration points of one integration segment and calculate
             # local mortar matrices
-            nsl = length(slave_element)
-            nm = length(master_element)
-            De = zeros(nsl, nsl)
-            Me = zeros(nsl, nm)
+            fill!(De, 0.0)
+            fill!(Me, 0.0)
             ge = zeros(field_dim*nsl)
             for ip in get_integration_points(slave_element, 2)
                 detJ = slave_element(ip, time, Val{:detJ})
@@ -162,20 +211,25 @@ function assemble!(problem::Problem{Mortar}, time::Real, ::Type{Val{1}})
                 xi = ip.coords[1]
                 xi_s = dot([1/2*(1-xi); 1/2*(1+xi)], xi1)
                 N1 = vec(get_basis(slave_element, xi_s, time))
+                Phi = Ae*N1
                 # project gauss point from slave element to master element in direction n_s
                 X_s = N1*X1 # coordinate in gauss point
                 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 
-                De += w*N1*N1'
-                Me += w*N1*N2'
+                De += w*Phi*N1'
+                Me += w*Phi*N2'
                 if props.adjust
+                    haskey(slave_element, "displacement") || continue
+                    haskey(master_element, "displacement") || continue
+                    norm(mean(X1) - X2[1]) / norm(X1[2] - X1[1]) < props.distval || continue
+                    norm(mean(X1) - X2[2]) / norm(X1[2] - X1[1]) < props.distval || continue
                     u1 = slave_element["displacement"](time)
                     u2 = master_element["displacement"](time)
                     x_s = X_s + N1*u1
                     x_m = X_m + N2*u2
-                    ge += w*vec((x_m-x_s)*N1')
+                    ge += w*vec((x_m-x_s)*Phi')
                 end
             end
 
@@ -199,6 +253,259 @@ function assemble!(problem::Problem{Mortar}, time::Real, ::Type{Val{1}})
 
 end
 
+# mesh tie 2d end
+
+# mesh tie 2d forwarddiff start
+
+function project_from_master_to_slave{E<:MortarElements2D}(
+    slave_element::Element{E}, x1_::DVTI, n1_::DVTI, x2::Vector, time::Float64;
+    tol=1.0e-10, max_iterations=20)
+
+    x1(xi1) = vec(get_basis(slave_element, [xi1], time))*x1_
+    dx1(xi1) = vec(get_dbasis(slave_element, [xi1], time))*x1_
+    n1(xi1) = vec(get_basis(slave_element, [xi1], time))*n1_
+    dn1(xi1) = vec(get_dbasis(slave_element, [xi1], time))*n1_
+    cross2(a, b) = cross([a; 0], [b; 0])[3]
+    R(xi1) = cross2(x1(xi1)-x2, n1(xi1))
+    dR(xi1) = cross2(dx1(xi1), n1(xi1)) + cross2(x1(xi1)-x2, dn1(xi1))
+
+    xi1 = 0.0
+    dxi1 = 0.0
+    for i=1:max_iterations
+        dxi1 = -R(xi1)/dR(xi1)
+        xi1 += dxi1
+        if norm(dxi1) < tol
+            return xi1
+        end
+    end
+
+    info("x1 = $(ForwardDiff.get_value(x1_.data))")
+    info("n1 = $(ForwardDiff.get_value(n1_.data))")
+    info("x2 = $(ForwardDiff.get_value(x2))")
+    info("xi1 = $(ForwardDiff.get_value(xi1)), dxi1 = $(ForwardDiff.get_value(dxi1))")
+    info("-R(xi1) = $(ForwardDiff.get_value(-R(xi1)))")
+    info("dR(xi1) = $(ForwardDiff.get_value(dR(xi1)))")
+    error("find projection from master to slave: did not converge")
+
+end
+
+function project_from_slave_to_master{E<:MortarElements2D}(
+    master_element::Element{E}, x1::Vector, n1::Vector, x2_::DVTI, time::Float64;
+    tol=1.0e-10, max_iterations=20)
+
+    x2(xi2) = vec(get_basis(master_element, [xi2], time))*x2_
+    dx2(xi2) = vec(get_dbasis(master_element, [xi2], time))*x2_
+    cross2(a, b) = cross([a; 0], [b; 0])[3]
+    R(xi2) = cross2(x2(xi2)-x1, n1)
+    dR(xi2) = cross2(dx2(xi2), n1)
+
+    xi2 = 0.0
+    dxi2 = 0.0
+    for i=1:max_iterations
+        dxi2 = -R(xi2) / dR(xi2)
+        xi2 += dxi2
+        if norm(dxi2) < tol
+            return xi2
+        end
+    end
+
+    error("find projection from slave to master: did not converge, last val: $xi2 and $dxi2")
+
+end
+
+""" 2d mesh tie using ForwardDiff.
+
+Construct .. + fc*la and C(d,la)=0
+
+"""
+function assemble!(problem::Problem{Mortar}, time::Float64, ::Type{Val{1}}, ::Type{Val{true}})
+
+    props = problem.properties
+    field_dim = get_unknown_field_dimension(problem)
+    field_name = get_parent_field_name(problem)
+    slave_elements = get_slave_elements(problem)
+    if field_name != "displacement"
+        error("mortar forwarddiff assembly: only displacement field with adjust=yes supported")
+    end
+
+    function calculate_interface(x::Vector)
+
+        ndofs = round(Int, length(x)/2)
+        nnodes = round(Int, ndofs/field_dim)
+        u = reshape(x[1:ndofs], field_dim, nnodes)
+        la = reshape(x[ndofs+1:end], field_dim, nnodes)
+        fc = zeros(u)
+        gap = zeros(u)
+        C = zeros(la)
+
+        S = Set{Int64}()
+        # 1. update nodal normals for slave elements
+        tangents = zeros(u)
+        for element in slave_elements
+            conn = get_connectivity(element)
+            push!(S, conn...)
+            X1 = element("geometry", time)
+            u1 = Field([u[:,i] for i in conn])
+            x1 = X1 + u1
+            dN = get_dbasis(element, [0.0], time)
+            tangent = sum([kron(dN[:,i], x1[i]') for i=1:length(x1)])
+            for nid in conn
+                tangents[:,nid] += tangent[:]
+            end
+        end
+
+        Q = [0.0 -1.0; 1.0 0.0]
+        normals = zeros(u)
+        for j in S
+            tangents[:,j] /= norm(tangents[:,j])
+            normals[:,j] = Q*tangents[:,j]
+        end
+
+        if props.rotate_normals
+            for j in S
+                normals[:,j] = -normals[:,j]
+            end
+        end
+
+        normals2 = Dict()
+        tangents2 = Dict()
+        for j in S
+            normals2[j] = normals[:,j]
+            tangents2[j] = tangents[:,j]
+        end
+        update!(slave_elements, "normal", time => normals2)
+        update!(slave_elements, "tangent", time => tangents2)
+
+        # 2. loop all slave elements
+        for slave_element in slave_elements
+
+            nsl = length(slave_element)
+            slave_element_nodes = get_connectivity(slave_element)
+            X1 = slave_element["geometry"](time)
+            u1 = Field(Vector[u[:,i] for i in slave_element_nodes])
+            x1 = X1 + u1
+            la1 = Field(Vector[la[:,i] for i in slave_element_nodes])
+            n1 = Field(Vector[normals[:,i] for i in slave_element_nodes])
+
+
+            # 3. loop all master elements
+            for master_element in slave_element["master elements"](time)
+
+                nm = length(master_element)
+                master_element_nodes = get_connectivity(master_element)
+                X2 = master_element["geometry"](time)
+                u2 = Field(Vector[u[:,i] for i in master_element_nodes])
+                x2 = X2 + u2
+
+                # 3.1 calculate segmentation
+                xi1a = project_from_master_to_slave(slave_element, x1, n1, x2[1], time)
+                xi1b = project_from_master_to_slave(slave_element, x1, n1, x2[2], time)
+#               xi1a = project_from_master_to_slave(slave_element, X2[1], time)
+#               xi1b = project_from_master_to_slave(slave_element, X2[2], time)
+                xi1 = clamp([xi1a; xi1b], -1.0, 1.0)
+                l = 1/2*abs(xi1[2]-xi1[1])
+                isapprox(l, 0.0) && continue # no contribution in this master element
+
+                # 3.2. bi-orthogonal basis
+                De = zeros(nsl, nsl)
+                Me = zeros(nsl, nsl)
+                Ae = zeros(nsl, nsl)
+                if props.dual_basis
+                    for ip in get_integration_points(slave_element, 3)
+                        detJ = slave_element(ip, time, Val{:detJ})
+                        w = ip.weight*detJ*l
+                        xi = ip.coords[1]
+                        xi_s = dot([1/2*(1-xi); 1/2*(1+xi)], xi1)
+                        N1 = vec(get_basis(slave_element, xi_s, time))
+                        De += w*diagm(N1)
+                        Me += w*N1*N1'
+                    end
+                    Ae = De*inv(Me)
+                else
+                    Ae = eye(nsl)
+                end
+
+                # 3.3. loop integration points of one integration segment and calculate
+                # local mortar matrices
+                for ip in get_integration_points(slave_element, 3)
+                    detJ = slave_element(ip, time, Val{:detJ})
+                    w = ip.weight*detJ*l
+                    #dN = get_dbasis(slave_element, ip, time)
+                    #j = sum([kron(dN[:,i], x1[i]') for i=1:length(x1)])
+                    #w = ip.weight*norm(j)*l
+                    
+                    xi = ip.coords[1]
+                    xi_s = dot([1/2*(1-xi); 1/2*(1+xi)], xi1)
+                    N1 = vec(get_basis(slave_element, xi_s, time))
+                    Phi = Ae*N1
+                    # project gauss point from slave element to master element in direction n_s
+                    x_s = N1*x1 # coordinate in gauss point
+                    n_s = N1*n1 # normal direction in gauss point
+                    #xi_m = project_from_slave_to_master(master_element, X_s, n_s, time)
+                    xi_m = project_from_slave_to_master(master_element, x_s, n_s, x2, time)
+                    N2 = vec(get_basis(master_element, xi_m, time))
+                    x_m = N2*x2 
+
+                    la_s = Phi*la1
+                    gn = dot(n_s, x_s-x_m)
+
+                    u_s = N1*u1
+                    u_m = N2*u2
+                    X_s = N1*X1
+                    X_m = N2*X2
+
+                    fc[:,slave_element_nodes] += w*la_s*N1'
+                    fc[:,master_element_nodes] -= w*la_s*N2'
+                    #gap[1,slave_element_nodes] += w*gn*Phi'
+                    gap[:,slave_element_nodes] += w*(u_s-u_m)*Phi'
+                    if props.adjust
+                        G = ForwardDiff.get_value(w*(X_s-X_m)*Phi')
+                        gap[:,slave_element_nodes] += G
+                    end
+                end
+
+            end # master elements done
+
+        end # slave elements done, contact virtual work ready
+
+        C = gap
+
+        info("interface residual ready")
+        return vec([fc C])
+
+    end
+
+    # x doesn't mean deformed configuration here
+    x = [problem.assembly.u; problem.assembly.la]
+    ndofs = round(Int, length(x)/2)
+    A, allresults = ForwardDiff.jacobian(calculate_interface, x,
+                            ForwardDiff.AllResults, cache=autodiffcache)
+    b = -ForwardDiff.value(allresults)
+
+    A = sparse(A)
+    b = sparse(b)
+    SparseMatrix.droptol!(A, 1.0e-12)
+    SparseMatrix.droptol!(b, 1.0e-12)
+
+    K = A[1:ndofs,1:ndofs]
+    C1 = transpose(A[1:ndofs,ndofs+1:end])
+    C2 = A[ndofs+1:end,1:ndofs]
+    D = A[ndofs+1:end,ndofs+1:end]
+    f = b[1:ndofs]
+    g = b[ndofs+1:end]
+
+    empty!(problem.assembly)
+    problem.assembly.K = K
+    problem.assembly.C1 = C1
+    problem.assembly.C2 = C2
+    problem.assembly.D = D
+    problem.assembly.f = f
+    problem.assembly.g = g
+
+end
+
+## 3d Mortar mesh tie
+
 function project_vertex_to_auxiliary_plane(p::Vector, x0::Vector, n0::Vector)
     return p - dot(p-x0, n0)*n0
 end
diff --git a/src/postprocess_utils.jl b/src/postprocess_utils.jl
index 26de9cd..ca33f40 100644
--- a/src/postprocess_utils.jl
+++ b/src/postprocess_utils.jl
@@ -1,42 +1,77 @@
 # This file is a part of JuliaFEM.
 # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
 
+using JuliaFEM
+
 """
 Calculate field values to nodal points from Gauss points using least-squares fitting.
 """
-function calc_nodal_values!(elements, field_name, field_dim, time)
-    A = SparseMatrixCOO()
-    b = SparseMatrixCOO()
-    for element in elements
-        gdofs = get_connectivity(element)
-        for ip in get_integration_points(element)
-            detJ = element(ip, time, Val{:detJ})
-            w = ip.weight*detJ
-            f = ip(field_name, time)
-            N = element(ip, time)
-            add!(A, gdofs, gdofs, w*kron(N', N))
-            for dim=1:field_dim
-                add!(b, gdofs, w*f[dim]*N, dim)
+function calc_nodal_values!(elements::Vector, field_name, field_dim, time;
+                            F=nothing, nz=nothing, b=nothing, return_F_and_nz=false)
+
+    if F == nothing
+        A = SparseMatrixCOO()
+        for element in elements
+            gdofs = get_connectivity(element)
+            for ip in get_integration_points(element)
+                detJ = element(ip, time, Val{:detJ})
+                w = ip.weight*detJ
+                N = element(ip, time)
+                add!(A, gdofs, gdofs, w*kron(N', N))
             end
         end
+        nz = get_nonzero_rows(A)
+        A = sparse(A)
+        A = 1/2*(A + A')
+        F = ldltfact(A[nz,nz])
     end
-    A = sparse(A)
-    b = sparse(b)
-    nz = get_nonzero_rows(A)
+
+    if b == nothing
+        b = SparseMatrixCOO()
+        for element in elements
+            gdofs = get_connectivity(element)
+            for ip in get_integration_points(element)
+                if !haskey(ip, field_name)
+                    info("warning: integration point does not have field $field_name")
+                    continue
+                end
+                detJ = element(ip, time, Val{:detJ})
+                w = ip.weight*detJ
+                f = ip(field_name, time)
+                N = element(ip, time)
+                for dim=1:field_dim
+                    add!(b, gdofs, w*f[dim]*N, dim)
+                end
+            end
+        end
+        b = sparse(b)
+    end
+
     x = zeros(size(b)...)
-    x[nz, :] = A[nz,nz] \ b[nz, :]
+    x[nz, :] = F \ b[nz, :]
     nodal_values = Dict()
     for i=1:size(x,1)
         nodal_values[i] = vec(x[i,:])
     end
-    update!(elements, field_name, nodal_values)
+    update!(elements, field_name, time => nodal_values)
+    if return_F_and_nz
+        return F, nz
+    end
+end
+
+function calc_nodal_values!(problem::Problem, field_name, field_dim, time)
+    # after all, it's just a mass matrix ...
+#   isempty(problem.assembly.M) && assemble!(problem, time, Val{:mass_matrix}; density=1.0, dual_basis=false, dim=1)
+#   M = sparse(problem.assembly.M)
+    # TODO: make test before implementation
+    calc_nodal_values!(problem.elements, field_name, field_dim, time)
 end
 
 """
 Return node ids + vector of values 
 """
 function get_nodal_vector(elements, field_name, time)
-    f = Dict{Int64, Vector{Float64}}()
+    f = Dict()
     for element in elements
         for (c, v) in zip(get_connectivity(element), element[field_name](time))
             if haskey(f, c)
@@ -50,3 +85,32 @@ function get_nodal_vector(elements, field_name, time)
     return node_ids, field
 end
 
+""" Update nodal field values from set of elements to another. Can be used to
+transform e.g. reaction force from boundary element set to surface of
+volume elements for easier postprocess.
+"""
+function copy_field!(src_elements::Vector, dst_elements::Vector, field_name, time)
+    dst_nodes = Set{Int64}()
+    for element in dst_elements
+        push!(dst_nodes, get_connectivity(element)...)
+    end
+    node_ids, field = get_nodal_vector(src_elements, field_name, time)
+    z = 0.0*first(field)
+    d = Dict()
+    for j in dst_nodes
+        d[j] = z
+    end
+    for (j, f) in zip(node_ids, field)
+        d[j] = f
+    end
+    for element in dst_elements
+        c = get_connectivity(element)
+        f = [d[j] for j in c]
+        update!(element, field_name, time => f)
+    end
+end
+
+function copy_field!(src_problem::Problem, dst_problem::Problem, field_name, time)
+    copy_field!(src_problem.elements, dst_problem.elements, field_name, time)
+end
+
diff --git a/src/postprocess_xdmf.jl b/src/postprocess_xdmf.jl
index 54f7159..8a6a09c 100644
--- a/src/postprocess_xdmf.jl
+++ b/src/postprocess_xdmf.jl
@@ -37,6 +37,9 @@ using JuliaFEM
 # > #define XDMF_3DRECTMESH     0x1101
 # > #define XDMF_3DCORECTMESH   0x1102
 
+get_xdmf_element_code(element::Element{Poi1})  = 0x0001
+get_xdmf_element_code(element::Element{Seg2})  = 0x0002
+get_xdmf_element_code(element::Element{Seg3})  = 0x0003
 get_xdmf_element_code(element::Element{Tri3})  = 0x0004
 get_xdmf_element_code(element::Element{Quad4}) = 0x0005
 get_xdmf_element_code(element::Element{Tet4})  = 0x0006
@@ -66,7 +69,7 @@ function XDMF()
     return XDMF(3, false, xdoc, domain, temporal_collection, Union{}, [])
 end
 
-function xdmf_new_result!(xdmf::XDMF, elements, time)
+function xdmf_new_result!(xdmf::XDMF, elements::Vector, time)
     grid = new_child(xdmf.temporal_collection, "Grid")
     set_attribute(grid, "Name", "Grid")
     time_ = new_child(grid, "Time")
@@ -128,9 +131,11 @@ function xdmf_new_result!(xdmf::XDMF, elements, time)
     add_text(dataitem, "\n"*join(s, "\n")*"\n")
 end
 
-function xdmf_save_field!(xdmf, elements, time, field_name; field_type="Scalar")
+function xdmf_save_field!(xdmf, elements::Vector, time, field_name; field_type="Scalar", debug=false)
     f = Dict()
+    field_dim = 0
     for element in elements
+        haskey(element, field_name) || continue
         g = element[field_name](time)
         conn = get_connectivity(element)
         for (i, c) in enumerate(conn)
@@ -139,10 +144,19 @@ function xdmf_save_field!(xdmf, elements, time, field_name; field_type="Scalar")
                 # paraview goes crazy if 2d model with 2d displacement vector
                 gi = [gi; 0.0]
             end
+            if field_dim == 0
+                field_dim = length(gi)
+            end
+            field_dim == length(gi) || error("several dimensions in field, dim = $field_dim.")
             f[c] = gi
         end
     end
 
+    if length(f) == 0
+        warn("xdmf_save_field!(): field $field_name was not found from set of elements")
+        return
+    end
+
     attribute = new_child(xdmf.current_grid, "Attribute")
     set_attribute(attribute, "Center", "Node")
     set_attribute(attribute, "Name", ucfirst(field_name))
@@ -151,16 +165,30 @@ function xdmf_save_field!(xdmf, elements, time, field_name; field_type="Scalar")
     set_attribute(dataitem, "DataType", "Float")
     set_attribute(dataitem, "Format", "XML")
     #set_attribute(dataitem, "Precision", 8)
+    debug && info("field dim = $field_dim")
+    debug && info(f)
     s = ASCIIString[]
     dim = 0
     for i in xdmf.permutation
-        push!(s, join(round(f[i], 5), " "))
-        dim += length(f[i])
+        gi = zeros(field_dim)
+        if haskey(f, i)
+            gi = f[i]
+        end
+        push!(s, join(round(gi, 5), " "))
+        dim += length(gi)
     end
     set_attribute(dataitem, "Dimensions", dim)
     add_text(dataitem, "\n"*join(s, "\n")*"\n")
 end
 
+function xdmf_save_field!(xdmf, problem::Problem, time, field_name; field_type="Scalar")
+    xdmf_save_field!(xdmf, problem.elements, time, field_name; field_type=field_type)
+end
+
+function xdmf_new_result!(xdmf, problem::Problem, time)
+    xdmf_new_result!(xdmf, problem.elements, time)
+end
+
 function xdmf_save!(xdmf, filename)
     save_file(xdmf.xdoc, filename)
 end
diff --git a/src/problems.jl b/src/problems.jl
index 465ed43..f1269d3 100644
--- a/src/problems.jl
+++ b/src/problems.jl
@@ -12,12 +12,15 @@ General linearized problem to solve
          C2*Δu +      D*λ = g
 """
 type Assembly
-    # for field assembly
+
     M :: SparseMatrixCOO  # mass matrix
-    K :: SparseMatrixCOO  # stiffness matrix
+
+    # for field assembly
+    K :: SparseMatrixCOO   # stiffness matrix
     Kg :: SparseMatrixCOO  # geometric stiffness matrix
-    f :: SparseMatrixCOO  # force vector
-#   f2 :: SparseMatrixCOO
+    f :: SparseMatrixCOO   # force vector
+    fg :: SparseMatrixCOO  # 
+
     # for boundary assembly
     C1 :: SparseMatrixCOO
     C2 :: SparseMatrixCOO
@@ -33,7 +36,6 @@ type Assembly
     la_prev :: Vector{Float64}  # previous solution vector u
     la_norm_change :: Real # change of norm in la
 
-    changed :: Bool  # flag to control is reassembly needed
 end
 
 function Assembly()
@@ -47,22 +49,38 @@ function Assembly()
         SparseMatrixCOO(),
         SparseMatrixCOO(),
         SparseMatrixCOO(),
+        SparseMatrixCOO(),
         [], [], Inf,
-        [], [], Inf,
-        true)
+        [], [], Inf)
 end
 
 function empty!(assembly::Assembly)
-    empty!(assembly.M)
     empty!(assembly.K)
     empty!(assembly.Kg)
     empty!(assembly.f)
+    empty!(assembly.fg)
     empty!(assembly.C1)
     empty!(assembly.C2)
     empty!(assembly.D)
     empty!(assembly.g)
     empty!(assembly.c)
-    assembly.changed = true
+end
+
+function isempty(assembly::Assembly)
+    T = isempty(assembly.K)
+    T &= isempty(assembly.Kg)
+    T &= isempty(assembly.f)
+    T &= isempty(assembly.fg)
+    T &= isempty(assembly.C1)
+    T &= isempty(assembly.C2)
+    T &= isempty(assembly.D)
+    T &= isempty(assembly.g)
+    T &= isempty(assembly.c)
+    return T
+end
+
+function get_dofs(assembly::Assembly)
+    return sort(unique(assembly.K.J))
 end
 
 type Problem{P<:AbstractProblem}
@@ -81,14 +99,15 @@ Examples
 --------
 Create vector-valued (dim=3) elasticity problem:
 
-julia> prob = Problem(Elasticity, "this is my problem", 3)
+julia> prob1 = Problem(Elasticity, "this is my problem", 3)
+julia> prob2 = Problem(Elasticity, 3)
 
 """
-function Problem{P<:FieldProblem}(::Type{P}, name::ASCIIString, dimension::Int64, elements=[], dofmap=Dict())
-    Problem{P}(name, dimension, "none", elements, dofmap, Assembly(), P())
+function Problem{P<:FieldProblem}(::Type{P}, name::ASCIIString, dimension::Int64)
+    Problem{P}(name, dimension, "none", [], Dict(), Assembly(), P())
 end
-function Problem{P<:FieldProblem}(::Type{P}, dimension::Int64, elements=[], dofmap=Dict())
-    Problem{P}("$P problem", dimension, "none", elements, dofmap, Assembly(), P())
+function Problem{P<:FieldProblem}(::Type{P}, dimension::Int64)
+    Problem{P}("$P problem", dimension, "none", [], Dict(), Assembly(), P())
 end
 
 """ Construct a new boundary problem.
@@ -100,14 +119,14 @@ Create Dirichlet boundary problem for vector-valued (dim=3) elasticity problem.
 julia> bc1 = Problem(Dirichlet, "support", 3, "displacement")
 
 """
-function Problem{P<:BoundaryProblem}(::Type{P}, name, dimension, parent_field_name, elements=[], dofmap=Dict())
-    Problem{P}(name, dimension, parent_field_name, elements, dofmap, Assembly(), P())
+function Problem{P<:BoundaryProblem}(::Type{P}, name, dimension, parent_field_name)
+    Problem{P}(name, dimension, parent_field_name, [], Dict(), Assembly(), P())
 end
-function Problem{P<:BoundaryProblem}(::Type{P}, main_problem::Problem, elements=[], dofmap=Dict())
+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)
-    Problem{P}(name, dimension, parent_field_name, elements, dofmap, Assembly(), P())
+    Problem{P}(name, dimension, parent_field_name, [], Dict(), Assembly(), P())
 end
 
 function get_formulation_type{P<:FieldProblem}(problem::Problem{P})
@@ -283,6 +302,14 @@ function get_gdofs(element::Element, dim::Int)
     return gdofs
 end
 
+function get_dofs(problem::Problem)
+    return get_dofs(problem.assembly)
+end
+
+function empty!(problem::Problem)
+    empty!(problem.assembly)
+end
+
 """ Return global degrees of freedom for element.
 
 Notes
diff --git a/src/solvers.jl b/src/solvers.jl
index 4b90d6d..5acd55f 100644
--- a/src/solvers.jl
+++ b/src/solvers.jl
@@ -7,49 +7,15 @@ type Solver{S<:AbstractSolver}
     name :: ASCIIString # some descriptive name for problem
     time :: Real        # current time
     problems :: Vector{Problem}
-    ndofs :: Int        # total dimension of global stiffness matrix, i.e., dim*nnodes
+    norms :: Vector{Tuple}  # solution norms for convergence studies
+    ndofs :: Int        # number of degrees of freedom in problem
     properties :: S
 end
 
-type Nonlinear <: AbstractSolver
-    iteration :: Int    # iteration counter
-    norms :: Vector{Tuple} # solution norms for convergence studies
-    min_iterations :: Int64
-    max_iterations :: Int64
-    convergence_tolerance :: Float64
-    error_if_no_convergence :: Bool
-    is_linear_system :: Bool # setting this to true makes assumption of one step convergence
-    linear_system_solver :: Symbol
-end
 
-function Nonlinear()
-    solver = Nonlinear(
-        0,     # iteration number
-        [],    # solution norms in (norm(u), norm(la)) tuples
-        1,     # min nonlinear iterations
-        10,    # max nonlinear iterations
-        5.0e-5, # nonlinear iteration convergence tolerance
-        true,  # throw error if no convergence
-        false, # is_linear_system
-        :DirectLinearSolver) # linear system solution method
-   return solver
-end
-
-function Solver{S<:AbstractSolver}(::Type{S}=Nonlinear,
-                                   name::ASCIIString="default solver",
-                                   time::Real=0.0, problems=[],
-                                   properties...)
+function Solver{S<:AbstractSolver}(::Type{S}, name="solver", properties...)
     variant = S(properties...)
-    solver = Solver{S}(name, time, problems, 0, variant)
-    return solver
-end
-
-""" For compatibility. """
-function Solver(name::ASCIIString="default solver",
-                time::Real=0.0, problems=[],
-                properties...)
-    variant = Nonlinear(properties...)
-    solver = Solver{Nonlinear}(name, time, problems, 0, variant)
+    solver = Solver{S}(name, 0.0, [], [], 0, variant)
     return solver
 end
 
@@ -72,50 +38,28 @@ end
 
 # one-liner helpers to identify problem types
 
-function is_field_problem(problem)
-    return false
-end
-function is_field_problem{P<:FieldProblem}(problem::Problem{P})
-    return true
-end
+is_field_problem(problem) = false
+is_field_problem{P<:FieldProblem}(problem::Problem{P}) = true
+is_boundary_problem(problem) = false
+is_boundary_problem{P<:BoundaryProblem}(problem::Problem{P}) = true
+get_field_problems(solver::Solver) = filter(is_field_problem, get_problems(solver))
+get_boundary_problems(solver::Solver) = filter(is_boundary_problem, get_problems(solver))
 
-function is_boundary_problem(problem)
-    return false
-end
-function is_boundary_problem{P<:BoundaryProblem}(problem::Problem{P})
-    return true
-end
+"""
+Posthook for field assembly. By default, do nothing.
+This can be used to make some modifications for assembly
+after all elements are assembled.
 
-function is_dirichlet_problem(problem)
-    return false
+Examples
+--------
+function field_assembly_posthook!(solver::Solver,
+                                  K::SparseMatrixCSC,
+                                  Kg::SparseMatrixCSC,
+                                  f::SparseMatrixCSC,
+                                  fg::SpareMatrixCSC)
+    info("doing stuff, size(K) = ", size(K))
 end
-function is_dirichlet_problem{P<:Problem{Dirichlet}}(problem::P)
-    return true
-end
-
-#=
-function is_mortar_problem{P<:Problem{Mortar}}(problem::P)
-    return true
-end
-=#
-
-function get_field_problems(solver::Solver)
-    filter(is_field_problem, solver.problems)
-end
-
-function get_boundary_problems(solver::Solver)
-    filter(is_boundary_problem, solver.problems)
-end
-
-function get_dirichlet_problems(solver::Solver)
-    filter(is_dirichlet_problem, solver.problems)
-end
-
-function get_mortar_problems(solver::Solver)
-    filter(is_mortar_problem, solver.problems)
-end
-
-""" Posthook for field assembly. By default, do nothing. """
+"""
 function field_assembly_posthook!
 end
 
@@ -127,7 +71,7 @@ solver :: Solver
 
 Returns
 -------
-K, f :: SparseMatrixCSC
+M, K, Kg, f, fg :: SparseMatrixCSC
 
 Notes
 -----
@@ -135,41 +79,41 @@ If several field problems exists, they are simply summed together, so
 problems must have unique node ids.
 
 """
-function get_field_assembly(solver::Solver; symmetric=true,
-                            with_mass_matrix=false,
-                            empty_after_append=true)
+function get_field_assembly(solver::Solver; show_info=true)
     problems = get_field_problems(solver)
+
     M = SparseMatrixCOO()
     K = SparseMatrixCOO()
     Kg = SparseMatrixCOO()
     f = SparseMatrixCOO()
+    fg = SparseMatrixCOO()
+
     for problem in problems
+        append!(M, problem.assembly.M)
         append!(K, problem.assembly.K)
         append!(Kg, problem.assembly.Kg)
         append!(f, problem.assembly.f)
-        with_mass_matrix && append!(M, problem.assembly.M)
-        empty_after_append && empty!(problem.assembly)
+        append!(fg, problem.assembly.fg)
     end
+
     if solver.ndofs == 0
         solver.ndofs = size(K, 1)
+        show_info && info("automatically determined problem dimension, ndofs = $(solver.ndofs)")
     end
+
+    M = sparse(M, solver.ndofs, solver.ndofs)
     K = sparse(K, solver.ndofs, solver.ndofs)
     Kg = sparse(Kg, solver.ndofs, solver.ndofs)
-    M = sparse(M, solver.ndofs, solver.ndofs)
-    if symmetric
-        K = 1/2*(K + K')
-        Kg = 1/2*(Kg + Kg')
-        M = 1/2*(M + M')
-    end
     f = sparse(f, solver.ndofs, 1)
+    fg = sparse(fg, solver.ndofs, 1)
 
     # run any posthook for assembly if defined
-    args = Tuple{Solver, SparseMatrixCSC, SparseMatrixCSC, SparseMatrixCSC}
+    args = Tuple{Solver, SparseMatrixCSC, SparseMatrixCSC, SparseMatrixCSC, SparseMatrixCSC}
     if method_exists(field_assembly_posthook!, args)
-        field_assembly_posthook!(solver, K, Kg, f)
+        field_assembly_posthook!(solver, K, Kg, fg, fg)
     end
 
-    return M, K, Kg, f
+    return M, K, Kg, f, fg
 end
 
 """ Posthook for boundary assembly. By default, do nothing. """
@@ -223,14 +167,13 @@ function get_boundary_assembly(solver::Solver)
         D += D_
         f += f_
         g += g_
-        empty!(problem.assembly)
     end
     return K, C1, C2, D, f, g
 end
 
 
 """
-Construct new basis such that u = P*uh + g
+Given C and g, construct new basis such that v = P*u + g
 
 Parameters
 ----------
@@ -262,7 +205,7 @@ 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}}; debug=false)
+function solve!(K, C1, C2, D, f, g, u, la, ::Type{Val{1}}; F=nothing, debug=false)
 
     nnz(D) == 0 || return false
     nz = get_nonzero_rows(C2)
@@ -289,57 +232,136 @@ function solve!(K, C1, C2, D, f, g, u, la, ::Type{Val{1}}; debug=false)
     end
 
     # solve interior domain using LDLt factorization
-    u[I] = ldltfact(K[I,I]) \ (f[I] - K[I,B]*u[B])
+    if F == nothing
+        F = ldltfact(K[I,I])
+    end
+    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 true
+    return F, 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}})
+function solve!(K, C1, C2, D, f, g, u, la, ::Type{Val{2}}; F=nothing)
     # 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))
-    x[nz] = lufact(A[nz,nz]) \ full(b[nz])
+    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 true
+    return F, true
 end
 
-function solve_linear_system(solver::Solver)
-    info("solving linear system of $(length(solver.problems)) problems.")
-    t0 = time()
+""" Default linear system solver for solver. """
+function solve_linear_system(solver::Solver; F=nothing, empty_assemblies_before_solution=true, show_info=true)
+    show_info && info("Solving problems ...")
+    t0 = Base.time()
 
-    # assemble field problems
-    M, K, Kg, f = get_field_assembly(solver)
-    # assemble boundary problems
+    # assemble field & boundary problems
+    # TODO: return same kind of set for both assembly types
+    # M1, K1, Kg1, f1, fg1, C11, C21, D1, g1 = get_field_assembly(solver)
+    # M2, K2, Kg2, f2, fg2, C12, C22, D2, g2 = get_boundary_assembly(solver)
+
+    M, K, Kg, f, fg = get_field_assembly(solver)
     Kb, C1, C2, D, fb, g = get_boundary_assembly(solver)
     K = K + Kg + Kb
+    f = f + fg + fb
     K = 1/2*(K + K')
-    f = f + fb
+    M = 1/2*(M + M')
+
+    # free up some memory before solution
+    for problem in get_problems(solver)
+        if empty_assemblies_before_solution
+            empty!(problem.assembly)
+        else
+            optimize!(problem.assembly)
+        end
+        gc()
+    end
 
     u = zeros(solver.ndofs)
     la = zeros(solver.ndofs)
 
     status = false
+    i = 0
     for i in [1, 2]
-        status = solve!(K, C1, C2, D, f, g, u, la, Val{i})
+        F, status = solve!(K, C1, C2, D, f, g, u, la, Val{i}; F=F)
         if status
-            info("succesfully solved Ax = b using solver #$i")
             break
         end
     end
     status || error("Failed to solve linear system!")
 
-    info("linear system solver: solved in ", time()-t0, " seconds. norm = ", norm(u))
-    return u, la
+    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
+end
+
+""" Default assembler for solver. """
+function assemble!(solver::Solver; show_info=true)
+    show_info && info("Assembling problems ...")
+    t0 = Base.time()
+    nproblems = 0
+    ndofs = 0
+    for problem in solver.problems
+        empty!(problem.assembly)
+        assemble!(problem, solver.time)
+        nproblems += 1
+        ndofs = max(ndofs, size(problem.assembly.K, 2))
+    end
+    solver.ndofs = ndofs
+    t1 = round(Base.time()-t0, 2)
+    show_info && info("Assembled $nproblems problems in $t1 seconds. ndofs = $ndofs.")
+end
+
+""" Default initializer for solver. """
+function initialize!(solver::Solver; show_info=true)
+    show_info && info("Initializing problems ...")
+    t0 = Base.time()
+    for problem in solver.problems
+        initialize!(problem, solver.time)
+    end
+    t1 = round(Base.time()-t0, 2)
+    show_info && info("Initialized problems in $t1 seconds.")
+end
+
+""" Default update for solver. """
+function update!(solver::Solver, u::Vector, la::Vector; show_info=true)
+    show_info && info("Updating problems ...")
+    t0 = Base.time()
+    for problem in solver.problems
+        u_new, la_new = update_assembly!(problem, u, la)
+        update_elements!(problem, u_new, la_new)
+    end
+    t1 = round(Base.time()-t0, 2)
+    show_info && info("Updated problems in $t1 seconds.")
+end
+
+### Nonlinear quasistatic solver
+
+type Nonlinear <: AbstractSolver
+    iteration :: Int        # iteration counter
+    min_iterations :: Int64 # minimum number of iterations
+    max_iterations :: Int64 # maximum number of iterations
+    convergence_tolerance :: Float64
+    error_if_no_convergence :: Bool # throw error if no convergence
+end
+
+function Nonlinear()
+    solver = Nonlinear(0, 1, 20, 5.0e-5, true)
+    return solver
 end
 
 """ Check convergence of problems.
@@ -348,7 +370,7 @@ Notes
 -----
 Default convergence criteria is obtained by checking each sub-problem convergence.
 """
-function has_converged(solver::Solver{Nonlinear};
+function has_converged(solver::Solver{Nonlinear}; show_info=false,
                        check_convergence_for_boundary_problems=false)
     properties = solver.properties
     converged = true
@@ -358,23 +380,24 @@ function has_converged(solver::Solver{Nonlinear};
         if is_field_problem(problem)
             has_converged = problem.assembly.u_norm_change < eps
             if isapprox(norm(problem.assembly.u), 0.0)
+                # trivial solution
                 has_converged = true
             end
-            info("Details for problem $(problem.name)")
-            info("Norm: $(norm(problem.assembly.u))")
-            info("Norm change: $(problem.assembly.u_norm_change)")
-            info("Has converged? $(has_converged)")
+            show_info && info("Details for problem $(problem.name)")
+            show_info && info("Norm: $(norm(problem.assembly.u))")
+            show_info && info("Norm change: $(problem.assembly.u_norm_change)")
+            show_info && info("Has converged? $(has_converged)")
         end
         if is_boundary_problem(problem) && check_convergence_for_boundary_problems
             has_converged = problem.assembly.la_norm_change/norm(problem.assembly.la) < eps
-            info("Details for problem $(problem.name)")
-            info("Norm: $(norm(problem.assembly.la))")
-            info("Norm change: $(problem.assembly.la_norm_change)")
-            info("Has converged? $(has_converged)")
+            show_info && info("Details for problem $(problem.name)")
+            show_info && info("Norm: $(norm(problem.assembly.la))")
+            show_info && info("Norm change: $(problem.assembly.la_norm_change)")
+            show_info && info("Has converged? $(has_converged)")
         end
         converged &= has_converged
     end
-    return converged || properties.is_linear_system
+    return converged
 end
 
 type NonlinearConvergenceError <: Exception
@@ -386,27 +409,8 @@ function Base.showerror(io::IO, exception::NonlinearConvergenceError)
     print(io, "nonlinear iteration did not converge in $max_iters iterations!")
 end
 
-function assemble!(solver::Solver; force_assembly=true)
-    info("Assembling problems ...")
-    tic()
-    for problem in solver.problems
-        if force_assembly # force reassembly
-            problem.assembly.changed = true
-        end
-        assemble!(problem, solver.time)
-    end
-    t1 = round(toq(), 2)
-    info("Assembled in $t1 seconds.")
-end
-
-function initialize!(solver::Solver)
-    for problem in solver.problems
-        initialize!(problem, solver.time)
-    end
-end
-
 """ Default solver for quasistatic nonlinear problems. """
-function call(solver::Solver{Nonlinear})
+function call(solver::Solver{Nonlinear}; show_info=true)
 
     properties = solver.properties
 
@@ -415,39 +419,105 @@ function call(solver::Solver{Nonlinear})
 
     # 2. start non-linear iterations
     for properties.iteration=1:properties.max_iterations
-        info("Starting nonlinear iteration #$(properties.iteration)")
+        show_info && info(repeat("-", 80))
+        show_info && info("Starting nonlinear iteration #$(properties.iteration)")
+        show_info && info("Increment time t=$(round(solver.time, 3))")
+        show_info && info(repeat("-", 80))
 
-        # 2.1 update linearized assemblies (if needed)
+        # 2.1 update linearized assemblies
         assemble!(solver)
 
-        # 2.2 call solver for linearized system (default: direct lu factorization)
-        info("Solve linear system ...")
-        tic()
-        u, la = solve_linear_system(solver)
-        push!(properties.norms, (norm(u), norm(la)))
-        t1 = round(toq(), 2)
-        info("Solved Ax = b in $t1 seconds.")
+        # 2.2 call solver for linearized system
+        F, u, la = solve_linear_system(solver)
 
         # 2.3 update solution back to elements
-        for problem in solver.problems
-            u_new, la_new = update_assembly!(problem, u, la)
-            update_elements!(problem, u_new, la_new)
-        end
+        update!(solver, u, la)
 
         # 2.4 check convergence
         if has_converged(solver)
             info("Converged in $(properties.iteration) iterations.")
-            if properties.iteration < properties.min_iterations
-                info("Converged but continuing")
-            else
-                return true
-            end
+            properties.iteration >= properties.min_iterations && return true
+            info("Convergence criteria met, but iteration < min_iterations, continuing...")
         end
     end
 
     # 3. did not converge
-    if properties.error_if_no_convergence
-        throw(NonlinearConvergenceError(solver))
+    properties.error_if_no_convergence && throw(NonlinearConvergenceError(solver))
+end
+
+""" Convenience function to call nonlinear solver. """
+function NonlinearSolver(problems...)
+    solver = Solver(Nonlinear, "default nonlinear solver")
+    if length(problems) != 0
+        push!(solver, problems...)
+    end
+    return solver
+end
+
+
+### Linear quasistatic solver
+
+""" Quasistatic solver for linear problems.
+
+Notes
+-----
+Main differences in this solver, compared to nonlinear solver are:
+1. system of problems is assumed to converge in one step
+2. reassembly of problem is done only if it's manually requested using empty!(problem.assembly)
+
+"""
+type Linear <: AbstractSolver
+    norms :: Vector{Tuple}
+end
+
+function Linear()
+    solver = Linear([])
+end
+
+function assemble!(solver::Solver{Linear}; show_info=true)
+    show_info && info("Assembling problems ...")
+    tic()
+    nproblems = 0
+    ndofs = 0
+    for problem in get_problems(solver)
+        if isempty(problem.assembly)
+            assemble!(problem, solver.time)
+            nproblems += 1
+        else
+            show_info && info("$(problem.name) already assembled, skipping.")
+        end
+        ndofs = max(ndofs, size(problem.assembly.K, 2))
+    end
+    solver.ndofs = ndofs
+    t1 = round(toq(), 2)
+    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)
+    t0 = Base.time()
+    show_info && info(repeat("-", 80))
+    show_info && info("Starting linear solver")
+    show_info && info("Increment time t=$(round(solver.time, 3))")
+    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)
+    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. """
+function LinearSolver(problems...)
+    solver = Solver(Linear, "default linear solver")
+    if length(problems) != 0
+        push!(solver, problems...)
+    end
+    return solver
+end
+
+### End of linear quasistatic solver
+
diff --git a/src/sparse.jl b/src/sparse.jl
index 2d09157..35a1bea 100644
--- a/src/sparse.jl
+++ b/src/sparse.jl
@@ -132,11 +132,12 @@ function add!(A::SparseMatrixCOO, dofs::Vector{Int}, data::Array{Float64}, dim::
     append!(A.V, vec(data))
 end
 
-""" Combine (I,J,V) values is possible. """
+""" Combine (I,J,V) values is possible to reduce memory usage. """
 function optimize!(A::SparseMatrixCOO)
     I, J, V = findnz(sparse(A))
-    A = SparseMatrixCOO(I, J, V)
-    gc()
+    A.I = I
+    A.J = J
+    A.V = V
 end
 
 """ Find all nonzero rows from sparse matrix.
@@ -163,6 +164,7 @@ function get_nonzero_columns(A::Union{SparseMatrixCOO, Matrix})
 end
 
 function size(A::SparseMatrixCOO)
+    isempty(A) && return (0, 0)
     return maximum(A.I), maximum(A.J)
 end
 
diff --git a/test/test_dirichlet.jl b/test/test_dirichlet.jl
index 59ddd42..900ffc8 100644
--- a/test/test_dirichlet.jl
+++ b/test/test_dirichlet.jl
@@ -31,6 +31,7 @@ Matrix([
     p1 = Problem(Dirichlet, "test problem 1", 1, "temperature")
     p1.properties.dual_basis = false
     p2 = Problem(Dirichlet, "test problem 2", 1, "temperature")
+    p2.properties.dual_basis = true
     assemble!(p1, element)
     assemble!(p2, element)
     C1 = full(p1.assembly.C1)
@@ -48,6 +49,7 @@ Matrix([
     p1 = Problem(Dirichlet, "quadratic 1", 1, "temperature")
     p1.properties.dual_basis = false
     p2 = Problem(Dirichlet, "quadratic 1", 1, "temperature")
+    p2.properties.dual_basis = true
     assemble!(p1, element)
     assemble!(p2, element)
     C1 = full(p1.assembly.C1)
@@ -112,3 +114,33 @@ end
 end
 =#
 
+@testset "test analytical boundary condition" begin
+    X = Dict{Int64, Vector{Float64}}(
+        1 => [0.0, 0.0],
+        2 => [1.0, 0.0])
+    element = Element(Seg2, [1, 2])
+    update!(element, "geometry", X)
+    update!(element, "displacement 1", 0.0)
+    f(xi, time) = begin
+        info("function call at xi = $xi, time = $time")
+        X = element("geometry", xi, time)
+        info("geometry at xi, X = $X")
+        val = X[1]*time
+        info("result for field at xi = $val")
+        return val
+    end
+    update!(element, "displacement 2", f)
+    p = Problem(Dirichlet, "test boundary", 2, "displacement")
+    push!(p, element)
+    assemble!(p, 0.0)
+    g1 = full(p.assembly.g, 4, 1)
+    @test isapprox(g1, [0.0, 0.0, 0.0, 0.0])
+    empty!(p.assembly)
+    assemble!(p, 1.0)
+    g2 = full(p.assembly.g, 4, 1)
+    C2 = full(p.assembly.C2, 4, 4)
+    u = C2 \ g2
+    info("u = $u")
+    @test isapprox(u, [0.0, 0.0, 0.0, 1.0])
+end
+
diff --git a/test/test_elasticity_2d_linear_with_surface_load.jl b/test/test_elasticity_2d_linear_with_surface_load.jl
index ddb0298..ee0648d 100644
--- a/test/test_elasticity_2d_linear_with_surface_load.jl
+++ b/test/test_elasticity_2d_linear_with_surface_load.jl
@@ -3,14 +3,17 @@
 
 using JuliaFEM
 using JuliaFEM.Preprocess
+using JuliaFEM.Postprocess
 using JuliaFEM.Test
+using JLD
 
-@testset "test 2d linear elasticity with surface load" begin
+function JuliaFEM.get_model(::Type{Val{Symbol("test 2d linear elasticity with surface + volume load")}})
     meshfile = "/geometry/2d_block/BLOCK_1elem.med"
     mesh = aster_read_mesh(Pkg.dir("JuliaFEM")*meshfile)
 
     # field problem
     block = Problem(Elasticity, "BLOCK", 2)
+    block.properties.store_fields = ["stress", "strain"]
     block.properties.formulation = :plane_stress
     block.properties.finite_strain = false
     block.properties.geometric_stiffness = false
@@ -34,6 +37,13 @@ using JuliaFEM.Test
 
     solver = Solver("solve block problem")
     push!(solver, block, bc_sym)
+    return solver
+end
+
+@testset "test 2d linear elasticity with surface + volume load" begin
+
+    solver = get_model("test 2d linear elasticity with surface + volume load")
+    block, bc_sym = solver.problems
     call(solver)
 
     f = 288.0
@@ -49,14 +59,39 @@ using JuliaFEM.Test
     for ip in get_integration_points(block.elements[1])
         eps = ip("strain")
         @printf "%i | %8.3f %8.3f | %8.3f %8.3f %8.3f\n" ip.id ip.coords[1] ip.coords[2] eps[1] eps[2] eps[3]
-        @test isapprox(eps, [u3; 0.0])
+        @test isapprox(eps, [u3[1], u3[2], 0.0])
     end
 
     info("stress")
     for ip in get_integration_points(block.elements[1])
         sig = ip("stress")
         @printf "%i | %8.3f %8.3f | %8.3f %8.3f %8.3f\n" ip.id ip.coords[1] ip.coords[2] sig[1] sig[2] sig[3]
-        @test isapprox(sig, [0.0; g; 0.0])
+        @test isapprox(sig, [0.0, g, 0.0])
     end
 
+    calc_nodal_values!(block.elements, "strain", 3, 0.0)
+    calc_nodal_values!(block.elements, "stress", 3, 0.0)
+    info(block.elements[1]["stress"](0.0))
+    node_ids, strain = get_nodal_vector(block.elements, "strain", 0.0)
+    node_ids, stress = get_nodal_vector(block.elements, "stress", 0.0)
+    @test isapprox(stress[1], [0.0, g, 0.0])
+    @test isapprox(strain[1], [u3[1], u3[2], 0.0])
 end
+
+@testset "test dump model to disk and read back before and after solution" begin
+    solver = get_model("test 2d linear elasticity with surface + volume load")
+    save("/tmp/model.jld", "linear_model", solver)
+    solver2 = load("/tmp/model.jld")["linear_model"]
+    call(solver2)
+    save("/tmp/model.jld", "results", solver2)
+    solver3 = load("/tmp/model.jld")["results"]
+    block = solver3["BLOCK"]
+    u3 = reshape(block.assembly.u, 2, 4)[:,3]
+    f = 288.0
+    g = 576.0
+    E = 288.0
+    nu = 1/3
+    u3_expected = f/E*[-nu, 1] + g/(2*E)*[-nu, 1]
+    @test isapprox(u3, u3_expected)
+end
+
diff --git a/test/test_elements.jl b/test/test_elements.jl
index d1a611f..170597a 100644
--- a/test/test_elements.jl
+++ b/test/test_elements.jl
@@ -59,7 +59,7 @@ end
 
 =#
 
-@testset "test add time dependent field to element" begin
+@testset "add time dependent field to element" begin
     el = Element(Seg2, [1, 2])
     u1 = Vector{Float64}[[0.0, 0.0], [0.0, 0.0]]
     u2 = Vector{Float64}[[1.0, 1.0], [1.0, 1.0]]
@@ -69,5 +69,34 @@ end
     @test isapprox(el("displacement", [0.0], 0.0), [0.0, 0.0])
     @test isapprox(el("displacement", [0.0], 0.5), [0.5, 0.5])
     @test isapprox(el("displacement", [0.0], 1.0), [1.0, 1.0])
+    el2 = Element(Poi1, [1])
+    update!(el2, "force 1", 0.0 => 1.0)
+end
+
+@testset "add CVTV field to element" begin
+    el = Element(Seg2, [1, 2])
+    f(xi, time) = xi[1]*time
+    update!(el, "my field", f)
+    v = el("my field", [1.0], 2.0)
+    @test isapprox(v, 2.0)
+end
+
+@testset "add DCTI to element" begin
+    el = Element(Quad4, [1, 2, 3, 4])
+    update!(el, "displacement load", DCTI([4.0, 8.0]))
+    @test isa(el["displacement load"], DCTI)
+    @test !isa(el["displacement load"].data, DCTI)
+    update!(el, "displacement load 2", [4.0, 8.0])
+    @test isa(el["displacement load 2"], DCTI)
+    update!(el, "temperature", [1.0, 2.0, 3.0, 4.0])
+    @test isa(el["temperature"], DVTI)
+end
+
+@testset "interpolate DCTI from element" begin
+    el = Element(Seg2, [1, 2])
+    update!(el, "foobar", 1.0)
+    fb = el("foobar", [0.0], 0.0)
+    @test isa(fb, Float64)
+    @test isapprox(fb, 1.0)
 end
 
diff --git a/test/test_fields.jl b/test/test_fields.jl
index 9ab554c..f33f059 100644
--- a/test/test_fields.jl
+++ b/test/test_fields.jl
@@ -25,3 +25,9 @@ end
     @test f.data == 2.0
 end
 
+@testset "test 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
diff --git a/test/test_mortar_2d_mesh_tie.jl b/test/test_mortar_2d_mesh_tie.jl
index 3c67c27..15dbbe3 100644
--- a/test/test_mortar_2d_mesh_tie.jl
+++ b/test/test_mortar_2d_mesh_tie.jl
@@ -45,7 +45,6 @@ end
     p1, p2, p3, p4 = get_test_model()
     p1.properties.formulation = :plane_stress
     p2.properties.formulation = :plane_stress
-    p4.properties.dimension = 1
     p4.properties.adjust = true
     p4.properties.rotate_normals = false
     solver = Solver(Nonlinear)
@@ -82,7 +81,6 @@ end
     interface_master_elements = create_elements(mesh, "UPPER_BOTTOM")
     update!(interface_slave_elements, "master elements", interface_master_elements)
     interface.elements = [interface_master_elements; interface_slave_elements]
-    interface.properties.dimension = 1
 
     solver = Solver()
     push!(solver, upper, lower, bc_upper, bc_lower, interface)
@@ -186,7 +184,7 @@ end
 
 function JuliaFEM.get_model(::Type{Val{Symbol("mesh tie with curved 2d block")}};
             dy=0.0, adjust=false, tolerance=0.0, rotate_normals=false, swap=false,
-            dual_basis=false)
+            dual_basis=false, use_forwarddiff=false)
 
     mesh = get_mesh("curved 2d block splitted to upper and lower")
 
@@ -221,9 +219,10 @@ function JuliaFEM.get_model(::Type{Val{Symbol("mesh tie with curved 2d block")}}
     update!(interface_slave_elements, "master elements", interface_master_elements)
     interface.elements = [interface_master_elements; interface_slave_elements]
     interface.properties.adjust = adjust
-    interface.properties.tolerance = tolerance
+    interface.properties.distval = tolerance
     interface.properties.rotate_normals = rotate_normals
     interface.properties.dual_basis = dual_basis
+    interface.properties.use_forwarddiff = use_forwarddiff
 
     solver = Solver(Nonlinear)
     push!(solver, upper, lower, bc_upper, bc_lower, interface)
@@ -276,4 +275,3 @@ end
     @test solver.properties.iteration == 2
     @test isapprox(norm(interface.assembly.u), 0.34318800698017704)
 end
-
diff --git a/test/test_mortar_2d_mesh_tie_forwarddiff.jl b/test/test_mortar_2d_mesh_tie_forwarddiff.jl
new file mode 100644
index 0000000..d2e7136
--- /dev/null
+++ b/test/test_mortar_2d_mesh_tie_forwarddiff.jl
@@ -0,0 +1,195 @@
+# 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.Preprocess
+using JuliaFEM.Postprocess
+using JuliaFEM.Test
+
+function JuliaFEM.get_mesh(::Type{Val{Symbol("curved 2d block splitted to upper and lower")}})
+    meshfile = Pkg.dir("JuliaFEM") * "/test/testdata/block_2d_curved.med"
+    mesh = aster_read_mesh(meshfile)
+end
+
+function JuliaFEM.get_model(::Type{Val{Symbol("mesh tie with curved 2d block")}};
+            dy=0.0, adjust=false, tolerance=0.0, rotate_normals=false, swap=false,
+            dual_basis=false, use_forwarddiff=true, finite_strain=false,
+            geometric_stiffness=false)
+
+    mesh = get_mesh("curved 2d block splitted to upper and lower")
+
+    upper = Problem(Elasticity, "upper", 2)
+    upper.properties.formulation = :plane_stress
+    upper.properties.finite_strain = finite_strain
+    upper.properties.geometric_stiffness = geometric_stiffness
+    upper.elements = create_elements(mesh, "UPPER")
+    update!(upper.elements, "youngs modulus", 96.0)
+    update!(upper.elements, "poissons ratio", 1/3)
+
+    lower = Problem(Elasticity, "lower", 2)
+    lower.properties.formulation = :plane_stress
+    lower.properties.finite_strain = finite_strain
+    lower.properties.geometric_stiffness = geometric_stiffness
+    lower.elements = create_elements(mesh, "LOWER")
+    update!(lower.elements, "youngs modulus", 96.0)
+    update!(lower.elements, "poissons ratio", 1/3)
+
+    bc_upper = Problem(Dirichlet, "upper boundary", 2, "displacement")
+    bc_upper.elements = create_elements(mesh, "UPPER_TOP")
+    update!(bc_upper.elements, "displacement 1", 0.0)
+    update!(bc_upper.elements, "displacement 2", dy)
+
+    bc_lower = Problem(Dirichlet, "lower boundary", 2, "displacement")
+    bc_lower.elements = create_elements(mesh, "LOWER_BOTTOM")
+    update!(bc_lower.elements, "displacement 1", 0.0)
+    update!(bc_lower.elements, "displacement 2", 0.0)
+
+    interface = Problem(Mortar, "interface between upper and lower block", 2, "displacement")
+    interface_slave_elements = create_elements(mesh, "LOWER_TOP")
+    interface_master_elements = create_elements(mesh, "UPPER_BOTTOM")
+    if swap
+        interface_slave_elements, interface_master_elements = interface_master_elements, interface_slave_elements
+    end
+    update!(interface_slave_elements, "master elements", interface_master_elements)
+    interface.elements = [interface_master_elements; interface_slave_elements]
+    interface.properties.adjust = adjust
+    interface.properties.distval = tolerance
+    interface.properties.rotate_normals = rotate_normals
+    interface.properties.dual_basis = dual_basis
+    interface.properties.use_forwarddiff = use_forwarddiff
+    interface.assembly.u = zeros(2*length(mesh.nodes))
+    interface.assembly.la = zeros(2*length(mesh.nodes))
+
+    solver = Solver(Nonlinear)
+    push!(solver, upper, lower, bc_upper, bc_lower, interface)
+
+    return solver
+
+end
+
+
+@testset "curved surface with adjust=true, standard lagrange, slave=lower surface, dy=0.0" begin
+    # TODO: analytical solution now known, verify using other fem software
+    solver = get_model("mesh tie with curved 2d block";
+        adjust=false, tolerance=10, dy=-0.1, rotate_normals=true,
+        dual_basis=true, use_forwarddiff=true, finite_strain=true,
+        geometric_stiffness=true)
+    call(solver)
+    interface = solver["interface between upper and lower block"]
+    @test solver.properties.iteration == 2
+    @test isapprox(norm(interface.assembly.u), 0.11339715157447851)
+end
+
+#=
+
+@testset "curved surface with adjust=true, dual lagrange, slave=lower surface, dy=0.0" begin
+    # TODO: analytical solution now known, verify using other fem software
+    solver = get_model("mesh tie with curved 2d block";
+        adjust=true, tolerance=10, dy=0.0, rotate_normals=true,
+        dual_basis=true, use_forwarddiff=true)
+    call(solver)
+    interface = solver["interface between upper and lower block"]
+    @test solver.properties.iteration == 2
+    # differs -- why?
+    @test isapprox(norm(interface.assembly.u), 0.11660422877751599)
+end
+
+@testset "curved surface with adjust=true, standard lagrange, slave=lower surface, dy=-0.1" begin
+    # TODO: analytical solution now known, verify using other fem software
+    solver = get_model("mesh tie with curved 2d block";
+        adjust=true, tolerance=10, dy=-0.1, rotate_normals=true,
+        dual_basis=false, use_forwarddiff=true)
+    call(solver)
+    interface = solver["interface between upper and lower block"]
+    @test solver.properties.iteration == 2
+    @test isapprox(norm(interface.assembly.u), 0.34230262165505887)
+end
+
+@testset "curved surface, adjust=true, dual basis, slave=lower surface, dy=-0.1" begin
+    # TODO: analytical solution now known, verify using other fem software
+    solver = get_model("mesh tie with curved 2d block";
+        adjust=true, tolerance=10, dy=-0.1, rotate_normals=true,
+        dual_basis=true, use_forwarddiff=true)
+    call(solver)
+    interface = solver["interface between upper and lower block"]
+    @test solver.properties.iteration == 2
+    @test isapprox(norm(interface.assembly.u), 0.34318800698017704)
+end
+
+=#
+
+function Base.isapprox(A::SparseMatrixCOO, B::SparseMatrixCOO)
+    A2 = sparse(A)
+    B2 = sparse(B, size(A2)...)
+    return isapprox(A2, B2)
+end
+
+function Base.isapprox(a1::Assembly, a2::Assembly)
+    T = isapprox(a1.K, a2.K)
+    T &= isapprox(a1.C1, a2.C1)
+    T &= isapprox(a1.C2, a2.C2)
+    T &= isapprox(a1.D, a2.D)
+    T &= isapprox(a1.f, a2.f)
+    T &= isapprox(a1.g, a2.g)
+    return T
+end
+
+@testset "compare forwarddiff solution to normal" begin
+    X = Dict(
+        1 => [0.0, 0.0],
+        2 => [1.0, 0.0],
+        3 => [0.0, 1.0],
+        4 => [1.0, 1.0])
+    u = Dict(
+        1 => [0.0, 0.0],
+        2 => [0.0, 0.0],
+        3 => [0.0, 0.0],
+        4 => [0.0, 0.0])
+    sel1 = Element(Seg2, [1, 2])
+    mel1 = Element(Seg2, [3, 4])
+    update!([sel1, mel1], "geometry", X)
+    update!([sel1, mel1], "displacement", u)
+    update!(sel1, "master elements", [mel1])
+    p1 = Problem(Mortar, "test 1", 2, "displacement")
+    p2 = Problem(Mortar, "test 2", 2, "displacement")
+    push!(p1, sel1, mel1)
+    push!(p2, sel1, mel1)
+    #p1.properties.adjust = true
+    p2.properties.use_forwarddiff = true
+    #p1.properties.dual_basis = true
+    #p2.properties.dual_basis = true
+    p2.assembly.u = zeros(8)
+    p2.assembly.la = zeros(8)
+    assemble!(p1, 0.0)
+    assemble!(p2, 0.0)
+    @test isapprox(p1.assembly, p2.assembly)
+
+    empty!(p1.assembly)
+    empty!(p2.assembly)
+    p1.properties.adjust = true
+    p2.properties.adjust = true
+    assemble!(p1, 0.0)
+    assemble!(p2, 0.0)
+    C11 = full(p1.assembly.C1, 4, 8)
+    C12 = full(p2.assembly.C1, 4, 8)
+    C21 = full(p1.assembly.C2, 4, 8)
+    C22 = full(p2.assembly.C2, 4, 8)
+    D1 = full(p1.assembly.D)
+    D2 = full(p2.assembly.D)
+    g1 = full(p1.assembly.g, 4, 1)
+    g2 = full(p2.assembly.g, 4, 1)
+    println("C1")
+    dump(C11)
+    dump(C12)
+    println("C2")
+    dump(C21)
+    dump(C22)
+    println("D")
+    dump(D1)
+    dump(D2)
+    println("g")
+    dump(g1)
+    dump(g2)
+    @test isapprox(p1.assembly, p2.assembly)
+end
+
diff --git a/test/test_xdmf.jl b/test/test_xdmf.jl
index 559ce30..8119e30 100644
--- a/test/test_xdmf.jl
+++ b/test/test_xdmf.jl
@@ -1,8 +1,6 @@
 # This file is a part of JuliaFEM.
 # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
 
-module XDMFTests
-
 using JuliaFEM
 using JuliaFEM.Postprocess
 using JuliaFEM.Test
@@ -126,4 +124,31 @@ function test_write_to_xml()
     end
 end
 
+@testset "write simple xmf file" 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])
+    u = Dict{Int64, Vector{Float64}}(
+        1 => [0.0, 0.0],
+        2 => [0.0, 0.0],
+        3 => [0.5, 1.0],
+        4 => [0.0, 0.0])
+    n = Dict{Int64, Vector{Float64}}(
+        2 => [1.0, 0.0],
+        3 => [1.0, 0.0])
+    el1 = Element(Quad4, [1, 2, 3, 4])
+    el2 = Element(Seg2, [2, 3])
+    update!([el1, el2], "geometry", X)
+    update!([el1, el2], "displacement", u)
+    update!(el2, "normal", n)
+    xdmf = XDMF()
+    xdmf.dimension = 2
+    xdmf_new_result!(xdmf, [el1, el2], 0.0)
+    xdmf_save_field!(xdmf, [el1, el2], 0.0, "displacement"; field_type="Vector")
+    xdmf_save_field!(xdmf, [el1, el2], 0.0, "normal"; field_type="Vector")
+    xdmf_save!(xdmf, "/tmp/test.xmf")
+    # TODO: how to test?
 end
+