diff --git a/src/JuliaFEM.jl b/src/JuliaFEM.jl
index d489f8b..0697ce0 100644
--- a/src/JuliaFEM.jl
+++ b/src/JuliaFEM.jl
@@ -25,7 +25,27 @@ export AbstractPoint, Point, IntegrationPoint, IP, Node
 include("elements.jl") # common element routines
 export Node, AbstractElement, Element, update!, get_connectivity, get_basis, get_dbasis
 include("lagrange_macro.jl") # Continuous Galerkin (Lagrange) elements generated using macro
-export Seg2, Seg3, Tri3, Tri6, Quad4, Hex8, Tet4, Tet10
+
+type Poi1 <: AbstractElement
+end
+
+function size(element::Element{Poi1})
+    return (0, 1)
+end
+
+function length(element::Element{Poi1})
+    return 1
+end
+
+function get_basis(element::Element{Poi1}, ip, time)
+    return [1]
+end
+
+function call(element::Element{Poi1}, ip, time, ::Type{Val{:detJ}})
+    return 1.0
+end
+
+export Poi1, Seg2, Seg3, Tri3, Tri6, Quad4, Hex8, Tet4, Tet10
 include("nurbs.jl")
 export NSeg, NSurf, NSolid, is_nurbs
 
@@ -80,10 +100,11 @@ export calculate_normals,
        calculate_normals!,
        project_from_slave_to_master,
        project_from_master_to_slave,
-       Mortar
+       Mortar, get_slave_elements
 
-### Contact mechanics ###
-#include("contact.jl")
+### Mortar methods, contact mechanics extension ###
+include("contact.jl")
+export Contact
 
 # rest of things
 include("utils.jl")
@@ -96,16 +117,31 @@ end
 
 module Preprocess
 include("preprocess.jl")
-export create_elements
+export create_elements, Mesh,
+       add_node!, add_nodes!,
+       add_element!, add_elements!,
+       add_element_to_element_set!,
+       add_node_to_node_set!,
+       find_nearest_nodes
 include("preprocess_abaqus_reader.jl")
 include("preprocess_abaqus_reader_old.jl")
 include("preprocess_aster_reader.jl")
 export aster_create_elements, parse_aster_med_file, is_aster_mail_keyword,
        parse_aster_header, aster_parse_nodes, aster_renumber_nodes!,
        aster_renumber_elements!, aster_combine_meshes, aster_read_mesh,
-       filter_by_element_set, filter_by_element_id
+       filter_by_element_set, filter_by_element_id, MEDFile
 end
 
+function get_mesh(mesh_name::ASCIIString, args...; kwargs...)
+    return get_mesh(Val{Symbol(mesh_name)}, args...; kwargs...)
+end
+
+function get_model(model_name::ASCIIString, args...; kwargs...)
+    return get_model(Val{Symbol(model_name)}, args...; kwargs...)
+end
+
+export get_mesh, get_model
+
 module Postprocess
 include("postprocess_utils.jl")
 export calc_nodal_values!, get_nodal_vector
@@ -120,6 +156,7 @@ if VERSION >= v"0.5-"
 else
     using BaseTestNext
 end
+
 export @test, @testset, @test_throws
 #include("test.jl")
 end
diff --git a/src/contact.jl b/src/contact.jl
index fa3d06f..7a23288 100644
--- a/src/contact.jl
+++ b/src/contact.jl
@@ -1,130 +1,208 @@
 # This file is a part of JuliaFEM.
 # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
 
-"""
-Currently two strategies exists:
-
-a) Remove inactive inequality constraints in element level. This is done in
-   assemble! if normal_condition is set to :Contact. For some reason this
-   leads to convergence issues.
-b) Remove inactive inequality constraints in assembly level. This is done in
-   posthook algorithm if inequality_constraints is set to true. This gives
-   more robust behavior.
-
-   Either use inequality_constraints=True OR :Contact + :Slip, but do not mix.
-
-   minimum_distance can be used to roughly skip integration of mortar
-   projections for elements that are "far enough" from each other. Increases
-   performance.
-
-"""
-type Mortar <: BoundaryProblem
-    formulation :: Symbol # :total, :incremental, :autodiff
-    dual_basis :: Bool
-    inequality_constraints :: Bool # Launch PDASS to solve inequality constraints
-    normal_condition :: Symbol # Tie or Contact
-    tangential_condition :: Symbol # Stick or Slip
-    maximum_distance :: Float64 # don't check for a contact if elements are far enough
-    store_debug_info :: Bool # for making debugging easier
-    always_inactive :: Vector{Int64}
-    always_in_contact :: Vector{Int64} # nodes in this list always in contact
-    always_in_stick :: Vector{Int64} # nodes in this list always in stick
-    always_in_slip :: Vector{Int64} # nodes in this list always in slip
-    contact :: Bool
-    friction :: Bool
-    gap_sign :: Int # gap sign convention
+type Contact <: BoundaryProblem
+    dimension :: Int
     rotate_normals :: Bool
+    finite_sliding :: Bool
+    friction :: Bool
+    dual_basis :: Bool
+    use_forwarddiff :: Bool
+    minimum_active_set_size :: Int
 end
 
-function Mortar()
-    Mortar(:total, true, false, :Tie, :Stick, Inf, false, [], [], [], [], false, false, -1, false)
+function Contact()
+    return Contact(-1, false, false, false, true, false, 0)
 end
 
-function get_unknown_field_name(::Type{Mortar})
+function get_unknown_field_name(problem::Problem{Contact})
     return "reaction force"
 end
 
-function get_formulation_type(problem::Problem{Mortar})
-    return problem.properties.formulation
+function get_formulation_type(problem::Problem{Contact})
+    return :incremental
 end
 
-macro debug(msg)
-    haskey(ENV, "DEBUG") || return
-    return msg
+typealias ContactElements2D Union{Seg2}
+
+function assemble!(problem::Problem{Contact}, time::Real)
+    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
+    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)
 end
 
-function assemble!(problem::Problem{Mortar}, time::Real)
-    elements = get_elements(problem)
-    if length(elements) == 0
-        info("$(typeof(problem)) : forget to add elements?")
-        return
+""" 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)
+
+    props = problem.properties
+    field_dim = get_unknown_field_dimension(problem)
+    field_name = get_parent_field_name(problem)
+    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)
+
+    # 2. loop all slave elements
+    for slave_element in slave_elements
+
+        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)
+        Q1_ = [n1[1] t1[1]]
+        Q2_ = [n1[2] t1[2]]
+        Z = zeros(2, 2)
+        Q2 = [Q1_ Z; Z Q2_]
+
+        # 3. loop all master elements
+        for master_element in slave_element["master elements"](time)
+
+            X2 = master_element["geometry"](time)
+            u2 = master_element["displacement"](time)
+            x2 = X2 + u2
+
+            # 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)
+            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'
+            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
+                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*Phi*N1'
+                Me += w*Phi*N2'
+                x_s = X_s + N1*u1
+                x_m = X_m + N2*u2
+                la_s = Phi*la1
+                ge += w*vec((x_m-x_s)*Phi')
+                lae += w*vec(la_s*Phi')
+            end
+
+            # add contribution to contact virtual work
+            sdofs = get_gdofs(problem, slave_element)
+            mdofs = get_gdofs(problem, master_element)
+            nsldofs = length(sdofs)
+            nmdofs = length(mdofs)
+            D2 = zeros(nsldofs, nsldofs)
+            M2 = zeros(nmdofs, nmdofs)
+            for i=1:field_dim
+                D2[i:field_dim:end, i:field_dim:end] += De
+                M2[i:field_dim:end, i:field_dim:end] += Me
+            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
+
+    end # slave elements done, contact virtual work ready
+
+    S = sort(collect(keys(normals))) # slave element nodes
+    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]
     end
-    # returns 3 if eldim 2 (tri3, quad4, ...) for 3d problems etc.
-    eldim = size(elements[1], 1)+1
-    assemble!(problem, time, Val{eldim})
-end
+    debug && info("Constitutive modelling ready")
 
-include("mortar_2d.jl")
-include("mortar_2d_autodiff.jl")
-include("mortar_3d.jl")
-include("mortar_3d_autodiff.jl")
-
-""" Remove inactive inequality constraints by using primal-dual active set strategy. """
-function boundary_assembly_posthook!(solver::Solver, problem::Problem{Mortar}, C1, C2, D, g)
-    problem.properties.inequality_constraints || return
-    info("PDASS: Starting primal-dual active set strategy to determine active constraints")
-    S = Set{Int64}()
-    for element in get_elements(problem)
-        haskey(element, "master elements") || continue
-        push!(S, get_connectivity(element)...)
-    end
-    S = sort(collect(S))
-    dim = get_unknown_field_dimension(problem)
-    ndofs = solver.ndofs
-    nnodes = round(Int, ndofs/dim)
-
-    c = reshape(full(problem.assembly.c, ndofs, 1), dim, nnodes)
-    A = find(c[1,:] .> 0)
-    A = intersect(A, S)
-    I = setdiff(S, A)
-
-    info("PDASS: contact nodes: $(sort(collect(S)))")
-    info("PDASS: active nodes: $(sort(collect(A)))")
-    info("PDASS: inactive nodes: $(sort(collect(I)))")
-
-    # remove any inactive nodes
-    for j in I
-        dofs = [dim*(j-1)+i for i=1:dim]
-        C1[dofs,:] = 0
-        C2[dofs,:] = 0
-        D[dofs,:] = 0
-        g[dofs,:] = 0
-    end
-
-    # handle tangential condition for active nodes
-    if problem.properties.tangential_condition == :Slip
-        for j in A
-            dofs = [dim*(j-1)+i for i=1:dim]
-            tangential_dofs = dofs[2:end]
-            D[tangential_dofs,dofs] = C2[tangential_dofs,dofs]
-            C2[tangential_dofs,:] = 0
-            g[tangential_dofs,:] = 0
+    # active / inactive node detection
+    A = Set()
+    I = Set()
+    la = problem.assembly.la
+    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")
+            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
+    problem.assembly.D = D
+    problem.assembly.g = g
 
     return
+
 end
 
-function assemble_prehook!(problem::Problem{Mortar}, time::Real)
-    info("mortar assemble prehook at time $time")
-    slaves = Set{Element}()
-    for element in get_elements(problem)
-        haskey(element, "master elements") || continue
-        push!(slaves, element)
-    end
-    info("$(length(slaves)) slave elements")
-    length(slaves) != 0 || error("no slave elements found for problem (forget to add masters?).")
-    info("mortar: update normal-tangential system.")
-    calculate_normal_tangential_coordinates!(collect(slaves), time)
-    info("mortar assemble prehook done.")
-end
diff --git a/src/elasticity.jl b/src/elasticity.jl
index ec8df38..9fbd435 100644
--- a/src/elasticity.jl
+++ b/src/elasticity.jl
@@ -157,11 +157,11 @@ function assemble{El<:Union{Tri3,Tri6,Quad4}}(problem::Problem{Elasticity}, elem
     return Km, Kg, f
 end
 
-function assemble{El<:Union{Seg2,Seg3}}(problem::Problem{Elasticity}, element::Element{El}, time::Real, ::Type{Val{:plane}})
+function assemble{El<:Union{Poi1,Seg2,Seg3}}(problem::Problem{Elasticity}, element::Element{El}, time::Real, ::Type{Val{:plane}})
 
     props = problem.properties
     dim = get_unknown_field_dimension(problem)
-    nnodes = size(element, 2)
+    nnodes = length(element)
     Km = zeros(dim*nnodes, dim*nnodes)
     Kg = zeros(dim*nnodes, dim*nnodes)
     f = zeros(dim*nnodes)
diff --git a/src/elements.jl b/src/elements.jl
index 734b832..49a04d8 100644
--- a/src/elements.jl
+++ b/src/elements.jl
@@ -67,7 +67,7 @@ function call(element::Element, field_name::ASCIIString, time)
     return element[field_name](time)
 end
 
-function call(element::Element, field_name::ASCIIString, ip, time)
+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)
diff --git a/src/integrate.jl b/src/integrate.jl
index e012a65..2aa8277 100644
--- a/src/integrate.jl
+++ b/src/integrate.jl
@@ -58,6 +58,12 @@ end
 
 ### "cartesian" elements, integration rules comes from tensor product
 
+### 0d elements
+
+function get_integration_points(element::Poi1)
+    [ (1.0, [] ) ]
+end
+
 ### 1d elements
 
 typealias CartesianLineElement Union{Seg2, Seg3, NSeg}
diff --git a/src/mortar.jl b/src/mortar.jl
index ddf97de..67b9cc8 100644
--- a/src/mortar.jl
+++ b/src/mortar.jl
@@ -6,13 +6,14 @@ type Mortar <: BoundaryProblem
     rotate_normals :: Bool
     adjust :: Bool
     tolerance :: Float64
+    dual_basis :: Bool
 end
 
 function Mortar()
-    return Mortar(-1, false, false, 0.0)
+    return Mortar(-1, false, false, 0.0, false)
 end
 
-function get_unknown_field_name(::Type{Mortar})
+function get_unknown_field_name(problem::Problem{Mortar})
     return "reaction force"
 end
 
@@ -38,7 +39,7 @@ function cross2(a, b)
     cross([a; 0], [b; 0])[3]
 end
 
-function get_slave_elements(problem::Problem{Mortar})
+function get_slave_elements(problem::Problem)
     filter(el -> haskey(el, "master elements"), get_elements(problem))
 end
 
@@ -109,8 +110,10 @@ function calculate_normals!(elements, time, ::Type{Val{1}}; rotate_normals=false
 end
 
 function assemble!(problem::Problem{Mortar}, time::Real)
-    if problem.dimension == -1
-        error("set interface dimension: problem.properties.dimension = 1 or 2")
+    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})
 end
diff --git a/src/postprocess_xdmf.jl b/src/postprocess_xdmf.jl
index a87a7eb..54f7159 100644
--- a/src/postprocess_xdmf.jl
+++ b/src/postprocess_xdmf.jl
@@ -134,7 +134,12 @@ function xdmf_save_field!(xdmf, elements, time, field_name; field_type="Scalar")
         g = element[field_name](time)
         conn = get_connectivity(element)
         for (i, c) in enumerate(conn)
-            f[c] = g[i]
+            gi = g[i]
+            if (field_type == "Vector") && (length(gi) < 3)
+                # paraview goes crazy if 2d model with 2d displacement vector
+                gi = [gi; 0.0]
+            end
+            f[c] = gi
         end
     end
 
diff --git a/src/preprocess.jl b/src/preprocess.jl
index 7d76948..aa40855 100644
--- a/src/preprocess.jl
+++ b/src/preprocess.jl
@@ -21,11 +21,30 @@ function add_node!(mesh::Mesh, nid::Int, ncoords::Vector{Float64})
     mesh.nodes[nid] = ncoords
 end
 
+function add_nodes!(mesh::Mesh, nodes::Dict{Int64, Vector{Float64}})
+    for (nid, ncoords) in nodes
+        add_node!(mesh, nid, ncoords)
+    end
+end
+
+function add_node_to_node_set!(mesh::Mesh, set_name::ASCIIString, nids...)
+    if !haskey(mesh.node_sets, set_name)
+        mesh.node_sets[set_name] = Set{Int64}()
+    end
+    push!(mesh.node_sets[set_name], nids...)
+end
+
 function add_element!(mesh::Mesh, elid::Int, eltype::Symbol, connectivity::Vector{Int64})
     mesh.elements[elid] = connectivity
     mesh.element_types[elid] = eltype
 end
 
+function add_elements!(mesh::Mesh, elements::Dict{Int64, Tuple{Symbol, Vector{Int64}}})
+    for (elid, (eltype, elcon)) in elements
+        add_element!(mesh, elid, eltype, elcon)
+    end
+end
+
 function add_element_to_element_set!(mesh::Mesh, set_name::ASCIIString, elids...)
     if !haskey(mesh.element_sets, set_name)
         mesh.element_sets[set_name] = Set{Int64}()
@@ -73,3 +92,15 @@ function create_elements(mesh::Mesh, element_set::ASCIIString)
     return create_elements(filter_by_element_set(mesh, element_set))
 end
 
+""" find npts nearest nodes form mesh and return id numbers as list. """
+function find_nearest_nodes(mesh::Mesh, coords::Vector, npts=1)
+    dist = Dict{Int64, Float64}()
+    for (nid, c) in mesh.nodes
+        dist[nid] = norm(coords-c)
+    end
+    s = sort(collect(dist), by=x->x[2])
+    nd = s[1:npts] # [(id1, dist1), (id2, dist2), ..., (id_npts, dist_npts)]
+    node_ids = [n[1] for n in nd]
+    return node_ids
+end
+
diff --git a/src/preprocess_aster_reader.jl b/src/preprocess_aster_reader.jl
index 31cce45..4555607 100644
--- a/src/preprocess_aster_reader.jl
+++ b/src/preprocess_aster_reader.jl
@@ -262,23 +262,38 @@ function get_mesh_names(med::MEDFile)
     return collect(keys(med.data["FAS"]))
 end
 
-function get_nodes(med::MEDFile, mesh_name)
+function get_nodes(med::MEDFile, nsets, mesh_name)
     increments = keys(med.data["ENS_MAA"][mesh_name])
     @assert length(increments) == 1
     increment = first(increments)
     nodes = med.data["ENS_MAA"][mesh_name][increment]["NOE"]
     node_ids = nodes["NUM"]
+    nset_ids = nodes["FAM"]
     nnodes = length(node_ids)
     node_coords = nodes["COO"]
     dim = round(Int, length(node_coords)/nnodes)
     node_coords = reshape(node_coords, nnodes, dim)'
-    d = Dict{Int64}{Vector{Float64}}()
+    d = Dict{Int64}{Tuple{Symbol, Vector{Float64}}}()
     for i=1:nnodes
-        d[node_ids[i]] = node_coords[:, i]
+        nset = Symbol(nsets[nset_ids[i]])
+        d[node_ids[i]] = (nset, node_coords[:, i])
     end
     return d
 end
 
+function get_node_sets(med::MEDFile, mesh_name)
+    ns = Dict{Int64, Symbol}(0 => :NALL)
+    haskey(med.data["FAS"][mesh_name], "NOEUD") || return ns
+    nsets = med.data["FAS"][mesh_name]["NOEUD"]
+    for nset in keys(nsets)
+        k = split(nset, "_")
+        nset_id = parse(Int, k[2])
+        nset_name = ascii(pointer(convert(Vector{UInt8}, nsets[nset]["GRO"]["NOM"][1])))
+        ns[nset_id] = Symbol(nset_name)
+    end
+    return ns
+end
+
 function get_element_sets(med::MEDFile, mesh_name)
     es = Dict{Int64, Symbol}()
     if !haskey(med.data["FAS"][mesh_name], "ELEME")
@@ -304,7 +319,8 @@ global const med_elmap = Dict{Symbol, Vector{Int}}(
     :QU4 => [1, 2, 3, 4],
     :HE8 => [4, 8, 7, 3, 1, 5, 6, 2],  # ..?
     :TE4 => [3, 2, 1, 4],
-    :T10 => [3, 2, 1, 4, 6, 5, 7, 10, 9, 8]
+    :T10 => [3, 2, 1, 4, 6, 5, 7, 10, 9, 8],
+    :PO1 => [1]
 #   :T10 => [3, 4, 1, 2, 10, 8, 7, 6, 9, 5]
 #   :T10 => [5, 9, 6, 7, 8, 10, 2, 1, 4, 3]
 )
@@ -354,7 +370,7 @@ Returns
 Dict containing fields "nodes" and "connectivity".
 
 """
-function parse_aster_med_file(fn::ASCIIString, mesh_name=nothing)
+function parse_aster_med_file(fn::ASCIIString, mesh_name=nothing; debug=false)
     med = MEDFile(fn)
     if isa(mesh_name, Void)
         mesh_names = get_mesh_names(med::MEDFile)
@@ -362,10 +378,15 @@ function parse_aster_med_file(fn::ASCIIString, mesh_name=nothing)
         length(mesh_names) == 1 || error("several meshes found from med, pick one: $all_meshes")
         mesh_name = mesh_names[1]
     end
+    nsets = get_node_sets(med, mesh_name)
     elsets = get_element_sets(med, mesh_name)
-    elset_names = join(values(elsets), ", ")
-    info("Found $(length(elsets)) element sets: $elset_names")
-    nodes = get_nodes(med, mesh_name)
+    if debug
+        elset_names = join(values(elsets), ", ")
+        info("Code Aster .med reader: found $(length(elsets)) element sets: $elset_names")
+        nset_names = join(values(nsets), ", ")
+        info("Code ASter .med reader: found $(length(nsets)) node sets: $nset_names")
+    end
+    nodes = get_nodes(med, nsets, mesh_name)
     conn = get_connectivity(med, elsets, mesh_name)
     result = Dict{ASCIIString, Any}()
     result["nodes"] = nodes
@@ -373,25 +394,33 @@ function parse_aster_med_file(fn::ASCIIString, mesh_name=nothing)
     return result
 end
 
+global const mapping = Dict(
+    :PO1 => :Poi1,
+    :SE2 => :Seg2,
+    :SE3 => :Seg3,
+    :TR3 => :Tri3,
+    :TR6 => :Tri6,
+    :QU4 => :Quad4,
+    :QU8 => :Quad8,
+    :QU9 => :Quad9,
+    :HE8 => :Hex8,
+    :H20 => :Hex20,
+    :TE4 => :Tet4,
+    :T10 => :Tet10)
+
 function aster_read_mesh(fn::ASCIIString, mesh_name=nothing)
     result = parse_aster_med_file(fn, mesh_name)
     mesh = Mesh()
-    for (nid, ncoords) in result["nodes"]
+    for (nid, (nset, ncoords)) in result["nodes"]
         add_node!(mesh, nid, ncoords)
+        add_node_to_node_set!(mesh, string(nset), nid)
     end
-    mapping = Dict(
-        :PO1 => :Poi1,
-        :SE2 => :Seg2,
-        :TR3 => :Tri3,
-        :TR6 => :Tri6,
-        :QU4 => :Quad4,
-        :HE8 => :Hex8,
-        :TE4 => :Tet4,
-        :T10 => :Tet10)
     for (elid, (eltype, elset, elcon)) in result["connectivity"]
+        haskey(mapping, eltype) || error("Code Aster .med reader: element type $eltype not found from mapping")
         add_element!(mesh, elid, mapping[eltype], elcon)
         add_element_to_element_set!(mesh, string(elset), elid)
     end
     return mesh
 end
 
+# TODO: refactor and remove obsolete stuff.
diff --git a/src/problems.jl b/src/problems.jl
index 8b9239b..465ed43 100644
--- a/src/problems.jl
+++ b/src/problems.jl
@@ -247,21 +247,14 @@ function update_elements!{P<:BoundaryProblem}(problem::Problem{P}, u, la)
     end
 end
 
-#=
-function add_postprocessor!(problem::Union{FieldProblem, BoundaryProblem}, postprocessor_name::Symbol, args...; kwargs...)
-    push!(problem.postprocessors, (postprocessor_name, args, kwargs))
-end
-
-function add_preprocessor!(problem::Union{FieldProblem, BoundaryProblem}, preprocessor_name::Symbol, args...; kwargs...)
-    push!(problem.preprocessors, (preprocessor_name, args, kwargs))
-end
-
-=#
-
 function get_elements(problem)
     return problem.elements
 end
 
+function update!(problem::Problem, field_name::ASCIIString, field)
+    update!(problem.elements, field_name, field)
+end
+
 """ Return the dimension of the unknown field of this problem. """
 function get_unknown_field_dimension(problem::Problem)
     return problem.dimension
diff --git a/src/solver_utils.jl b/src/solver_utils.jl
index 4987fa8..48c4197 100644
--- a/src/solver_utils.jl
+++ b/src/solver_utils.jl
@@ -335,9 +335,11 @@ function handle_overconstraint_error!(problem, nodes, all_dofs, C1_, C1, C2_, C2
             continue
         end
            
-        show_info && info("unable to resolve overconstrained situation, not continuing")
+        info("unable to resolve overconstrained situation, not continuing")
+        show_rows_in_constraint_matrix(dofs, C2, D; show_status=false)
+        show_rows_in_constraint_matrix(dofs, C2_, D_; show_status=false)
+        show_related_equations(dofs, C2, C2_, D, D_)
         throw("failed to resolve overconstraint situation")
-        show_info && info()
     end
 end
 
diff --git a/src/solvers.jl b/src/solvers.jl
index 6bac12a..4b90d6d 100644
--- a/src/solvers.jl
+++ b/src/solvers.jl
@@ -53,10 +53,23 @@ function Solver(name::ASCIIString="default solver",
     return solver
 end
 
+function get_problems(solver::Solver)
+    return solver.problems
+end
+
 function push!(solver::Solver, problem)
     push!(solver.problems, problem)
 end
 
+function getindex(solver::Solver, problem_name::ASCIIString)
+    for problem in get_problems(solver)
+        if problem.name == problem_name
+            return problem
+        end
+    end
+    throw(KeyError(problem_name))
+end
+
 # one-liner helpers to identify problem types
 
 function is_field_problem(problem)
@@ -231,8 +244,11 @@ function create_projection(C::SparseMatrixCSC, g; S=nothing, tol=1.0e-12)
     end
     # FIXME: this creates dense matrices
     # efficiency / memory usage is a question
-    P = sparse(C[S,:] \ full(C[S,:]))
-    h = sparse(C[S,:] \ full(g[S]))
+    M = get_nonzero_columns(C)
+    F = qrfact(C[S,:])
+    P = spzeros(n,m)
+    P[:,M] = sparse(F \ full(C[S,M]))
+    h = sparse(F \ full(g[S]))
     resize!(P, n, m)
     resize!(h, n, 1)
     P = speye(n) - P
diff --git a/test/runtests.jl b/test/runtests.jl
index 1516239..9f34548 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,19 +1,23 @@
 # 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.Test
 
 function run_tests(; quiet=false)
 
-    test_files = readdir(Pkg.dir("JuliaFEM")*"/test")
-    test_files = filter(f -> (startswith(f, "test_") & endswith(f, ".jl")), test_files)
-    for test_file in test_files
-        if !quiet
-            info("Running tests from file $test_file")
+    maybe_test_files = readdir(Pkg.dir("JuliaFEM")*"/test")
+    is_test_file(fn) = startswith(fn, "test_") & endswith(fn, ".jl")
+    test_files = filter(is_test_file, maybe_test_files)
+    #test_files = ["test_nodal_constraints.jl"]
+
+    body = quote
+        @testset "JuliaFEM" begin
+            for fn in $test_files
+                @testset "$fn" begin include(fn) end
+            end
         end
-        include(test_file)
     end
+    eval(body)
 
 end
 
diff --git a/test/test_contact_2d_small_sliding.jl b/test/test_contact_2d_small_sliding.jl
new file mode 100644
index 0000000..a78b8b4
--- /dev/null
+++ b/test/test_contact_2d_small_sliding.jl
@@ -0,0 +1,147 @@
+# 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
+import JuliaFEM: get_mesh, get_model
+
+function get_mesh(::Type{Val{Symbol("curved 2d mesh model")}})
+    meshfile = Pkg.dir("JuliaFEM") * "/test/testdata/block_2d_curved.med"
+    mesh = aster_read_mesh(meshfile)
+end
+
+function get_model(::Type{Val{Symbol("curved 2d contact small sliding")}})
+    
+    mesh = get_mesh("curved 2d mesh model")
+
+    upper = Problem(Elasticity, "upper", 2)
+    upper.properties.formulation = :plane_stress
+    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.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", -0.15)
+
+    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(Contact, "contact between upper and lower block", 2, "displacement")
+    interface.properties.dimension = 1
+    interface.properties.rotate_normals = true
+    interface_slave_elements = create_elements(mesh, "LOWER_TOP")
+    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]
+
+    solver = Solver(Nonlinear)
+    push!(solver, upper, lower, bc_upper, bc_lower, interface)
+    return solver
+
+end
+
+@testset "test all nodes in contact" begin
+    # FIXME: needs verification of some other fem software
+    solver = get_model("curved 2d contact small sliding")
+    call(solver)
+    upper, lower, bc_upper, bc_lower, interface = solver.problems
+    @test isapprox(norm(interface.assembly.u), 0.49563347601324315)
+end
+
+
+function get_mesh(::Type{Val{Symbol("hertz contact, full 2d model")}})
+    meshfile = Pkg.dir("JuliaFEM") * "/test/testdata/hertz_2d_full.med"
+    mesh = aster_read_mesh(meshfile)
+end
+
+function get_model(::Type{Val{Symbol("hertz contact, full 2d model")}})
+    # from fenet d3613 advanced finite element contact benchmarks
+    # a = 6.21 mm, pmax = 3585 MPa
+    # this is a very dense mesh and for that reason pmax is not very
+    # (only 6 elements in -20 .. 20 mm contact zone, 3 elements in contact
+    # instead integrate pressure in normal and tangential direction
+    mesh = get_mesh("hertz contact, full 2d model")
+
+    upper = Problem(Elasticity, "CYLINDER", 2)
+    upper.properties.formulation = :plane_strain
+    upper.elements = create_elements(mesh, "CYLINDER")
+    update!(upper.elements, "youngs modulus", 70.0e3)
+    update!(upper.elements, "poissons ratio", 0.3)
+
+    lower = Problem(Elasticity, "BLOCK", 2)
+    lower.properties.formulation = :plane_strain
+    lower.elements = create_elements(mesh, "BLOCK")
+    update!(lower.elements, "youngs modulus", 210.0e3)
+    update!(lower.elements, "poissons ratio", 0.3)
+
+    # support block to ground
+    bc_fixed = Problem(Dirichlet, "fixed", 2, "displacement")
+    bc_fixed.elements = create_elements(mesh, "FIXED")
+    update!(bc_fixed.elements, "displacement 2", 0.0)
+
+    # symmetry line
+    bc_sym_23 = Problem(Dirichlet, "symmetry line 23", 2, "displacement")
+    bc_sym_23.elements = create_elements(mesh, "SYM23")
+    update!(bc_sym_23.elements, "displacement 1", 0.0)
+
+    nid = find_nearest_nodes(mesh, [0.0, 100.0])
+    #load = Problem(Dirichlet, "load", 2, "displacement")
+    load = Problem(Elasticity, "point load", 2)
+    load.properties.formulation = :plane_strain
+    load.elements = [Element(Poi1, nid)]
+    #update!(load.elements, "displacement 2", -10.0)
+    update!(load.elements, "displacement traction force 2", -35.0e3)
+
+    contact = Problem(Contact, "contact between block and cylinder", 2, "displacement")
+    contact.properties.rotate_normals = true
+    contact_slave_elements = create_elements(mesh, "CYLINDER_TO_BLOCK")
+    contact_master_elements = create_elements(mesh, "BLOCK_TO_CYLINDER")
+    update!(contact_slave_elements, "master elements", contact_master_elements)
+    contact.elements = [contact_master_elements; contact_slave_elements]
+    
+    solver = Solver(Nonlinear)
+    push!(solver, upper, lower, bc_fixed, bc_sym_23, load, contact)
+    return solver
+
+end
+
+@testset "test frictionless hertz contact, 2d plane strain" begin
+    solver = get_model("hertz contact, full 2d model")
+    call(solver)
+    upper, lower, bc_fixed, bc_sym_23, load, contact = solver.problems
+    slaves = get_slave_elements(contact)
+    node_ids, la = get_nodal_vector(slaves, "reaction force", 0.0)
+    node_ids, n = get_nodal_vector(slaves, "normal", 0.0)
+    pres = [dot(ni, lai) for (ni, lai) in zip(n, la)]
+    @test isapprox(maximum(pres), 4060.010799583303)
+    # integrate pressure in normal and tangential direction
+    Rn = 0.0
+    Rt = 0.0
+    Q = [0.0 -1.0; 1.0 0.0]
+    time = 0.0
+    for sel in slaves
+        for ip in get_integration_points(sel)
+            w = ip.weight*sel(ip, time, Val{:detJ})
+            n = sel("normal", ip, time)
+            t = Q'*n
+            la = sel("reaction force", ip, time)
+            Rn += w*dot(n, la)
+            Rt += w*dot(t, la)
+        end
+    end
+    @test isapprox(Rn, 35.0e3; rtol=0.0015)
+    @test isapprox(Rt, 0.0; atol=10.0)
+end
+
diff --git a/test/test_mortar_2d_contact.jl b/test/test_mortar_2d_contact.jl
index e056cfb..4cf86f3 100644
--- a/test/test_mortar_2d_contact.jl
+++ b/test/test_mortar_2d_contact.jl
@@ -1,263 +1,85 @@
 # 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.Test
 
-using JuliaFEM.Core: Node, Seg2, Tri3, update!, calculate_normal_tangential_coordinates!,
-                     PlaneStressLinearElasticityProblem, DirichletProblem, MortarProblem,
-                     get_elements, DirectSolver, calculate_nodal_vector, FieldAssembly,
-                     FieldProblem, set_linear_system_solver!, set_nonlinear_max_iterations!,
-                     get_elements, Element, BoundaryAssembly, BoundaryProblem,
-                     get_integration_points, get_jacobian, get_connectivity, StandardBasis,
-                     add_postprocessor!, add_preprocessor!, SparseMatrixCOO, add!,
-                     add_linear_system_solver_preprocessor!,
-                     add_linear_system_solver_postprocessor!
-
-import JuliaFEM.Core: assemble_preprocess!, assemble_postprocess!,
-                      linear_system_solver_preprocess!, linear_system_solver_postprocess!
-
-macro debug(msg)
-    haskey(ENV, "DEBUG") || return
-    return msg
+function JuliaFEM.get_mesh(::Type{Val{Symbol("two elements 1.0x0.5 with 0.1 gap in y direction")}})
+    mesh = Mesh()
+    add_node!(mesh, 1, [0.0, 0.0])
+    add_node!(mesh, 2, [1.0, 0.0])
+    add_node!(mesh, 3, [1.0, 0.5])
+    add_node!(mesh, 4, [0.0, 0.5])
+    add_node!(mesh, 5, [0.0, 0.6])
+    add_node!(mesh, 6, [1.0, 0.6])
+    add_node!(mesh, 7, [1.0, 1.1])
+    add_node!(mesh, 8, [0.0, 1.1])
+    add_element!(mesh, 1, :Quad4, [1, 2, 3, 4])
+    add_element!(mesh, 2, :Quad4, [5, 6, 7, 8])
+    add_element!(mesh, 3, :Seg2, [1, 2])
+    add_element!(mesh, 4, :Seg2, [7, 8])
+    add_element!(mesh, 5, :Seg2, [4, 3])
+    add_element!(mesh, 6, :Seg2, [6, 5])
+    add_element_to_element_set!(mesh, "LOWER", 1)
+    add_element_to_element_set!(mesh, "UPPER", 2)
+    add_element_to_element_set!(mesh, "LOWER_BOTTOM", 3)
+    add_element_to_element_set!(mesh, "UPPER_TOP", 4)
+    add_element_to_element_set!(mesh, "LOWER_TOP", 5)
+    add_element_to_element_set!(mesh, "UPPER_BOTTOM", 6)
+    return mesh
 end
 
-function calculate_normal_tangential_coordinates(elements::Vector{Element}, time::Real)
-    P = SparseMatrixCOO()
-    field_dim = 2
-    for element in elements
-        haskey(element, "normal-tangential coordinates") || continue
-        for ip in get_integration_points(element, Val{2})
-            J = get_jacobian(element, ip, time)
-            w = ip.weight*norm(J)
-            nt = transpose(element("normal-tangential coordinates", ip, time))
-            normal = nt[1,:]
-            tangent = nt[2,:]
-            for nid in get_connectivity(element)
-                ndofs = [2*(nid-1)+1, 2*(nid-1)+2]
-                add!(P, [2*(nid-1)+1], ndofs, normal)
-                add!(P, [2*(nid-1)+2], ndofs, tangent)
-            end
-        end
-    end
-    P = sparse(P)
-    for i=1:size(P,1)
-        n = norm(P[i,:])
-        if n > 0.0
-            P[i,:] = P[i,:] / n
-        end
-    end
-    return SparseMatrixCOO(P)
+function JuliaFEM.get_model(::Type{Val{Symbol("two element contact")}})
+
+    mesh = get_mesh("two elements 1.0x0.5 with 0.1 gap in y direction")
+
+    upper = Problem(Elasticity, "UPPER", 2)
+    upper.properties.formulation = :plane_stress
+    upper.elements = create_elements(mesh, "UPPER")
+    update!(upper.elements, "youngs modulus", 288.0)
+    update!(upper.elements, "poissons ratio", 1/3)
+
+    lower = Problem(Elasticity, "LOWER", 2)
+    lower.properties.formulation = :plane_stress
+    lower.elements = create_elements(mesh, "LOWER")
+    update!(lower.elements, "youngs modulus", 288.0)
+    update!(lower.elements, "poissons ratio", 1/3)
+
+    bc_upper = Problem(Dirichlet, "UPPER_TOP", 2, "displacement")
+    bc_upper.elements = create_elements(mesh, "UPPER_TOP")
+    #update!(bc_upper.elements, "displacement 1", -17/90)
+    #update!(bc_upper.elements, "displacement 1", -17/90)
+    update!(bc_upper.elements, "displacement 1", -0.2)
+    update!(bc_upper.elements, "displacement 2", -0.2)
+
+    bc_lower = Problem(Dirichlet, "LOWER_BOTTOM", 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(Contact, "LOWER_TO_UPPER", 2, "displacement")
+    interface.properties.dimension = 1
+    interface_slave_elements = create_elements(mesh, "LOWER_TOP")
+    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]
+
+    solver = Solver(Nonlinear)
+    push!(solver, upper, lower, bc_upper, bc_lower, interface)
+    return solver
+
 end
 
-function assemble_postprocess!(assembly, problem, time::Real, ::Type{Val{:remove_constraint_from_dofs}}, dofs)
-    # giving additional arguments and keywords is possible too
-    C1 = sparse(assembly.C1)
-    C2 = sparse(assembly.C2)
-    info("removing constraints from dofs: $dofs")
-    for d in dofs
-        C1[d,:] = 0
-        C2[d,:] = 0
-    end
-    assembly.C1 = C1
-    assembly.C2 = C2
+@testset "test simple two element contact" begin
+    solver = get_model("two element contact")
+    call(solver)
+    contact = solver["LOWER_TO_UPPER"]
+    master = first(contact.elements)
+    slave = last(contact.elements)
+    u = master("displacement", [0.0], 0.0)
+    la = slave("reaction force", [0.0], 0.0)
+    info("u = $u, la = $la")
+    @test isapprox(u, [-0.2, -0.15])
+    @test isapprox(la, [0.0, 30.375])
 end
-
-function assemble_postprocess!(assembly, problem, time::Real, ::Type{Val{:remove_tangential_constraints}})
-    # example how to use postprocessor to manipulate constraint matrix before summing assemblies together
-    info("postprocess mortar assembly: remove contraints in tangent direction on boundary.")
-    dim = 12
-    C1 = sparse(assembly.C1, dim, dim)
-    C2 = sparse(assembly.C2, dim, dim)
-    P = calculate_normal_tangential_coordinates(get_elements(problem), time)
-    P = sparse(P, dim, dim)
-    info("projection matrix for normals: ")
-    dump(round(full(P), 3))
-    C1 = P*C1
-    C2 = P*C2
-    for i=2:2:dim
-        C1[i,:] = 0
-        C2[i,:] = 0
-    end
-    assembly.C1 = C1
-    assembly.C2 = C2
-    info("postprocess mortar assembly: done.")
-end
-
-function assemble_postprocess!(assembly, problem, time::Real, ::Type{Val{:primal_dual_active_set_strategy}})
-    info("PDASS: determining active contact set")
-    dim = 12
-    C1 = sparse(assembly.C1, dim, dim)
-    C2 = sparse(assembly.C2, dim, dim)
-
-    info("PDASS: constraint matrix C1")
-#   dump(round(full(C1[1:2:end,1:2:end]), 3))
-    dump(round(full(C1), 3))
-    info("PDASS: constraint matrix C2")
-#   dump(round(full(C2[1:2:end,1:2:end]), 3))
-    dump(round(full(C2), 3))
-
-    elements = get_elements(problem)
-    P = calculate_normal_tangential_coordinates(get_elements(problem), time)
-    P = sparse(P, dim, dim)
-
-    la = calculate_nodal_vector("reaction force", 2, elements, time)
-    X = calculate_nodal_vector("geometry", 2, elements, time)
-    u = calculate_nodal_vector("displacement", 2, elements, time)
-    resize!(la, 12)
-    resize!(X, 12)
-    resize!(u, 12)
-    x = X+u
-    info("x")
-    dump(reshape(round(x, 2), 2, 6))
-    P = sparse(eye(dim))
-    C1 = P*C1
-    C2 = P*C2
-    gn = P*C1*X #*4/6 ..?
-    un = P*C1*u
-    la = P*la
-    info("weighted gap in nt =")
-    dump(reshape(round(gn, 2), 2, 6))
-    info("weighted u in nt =")
-    dump(reshape(round(un, 2), 2, 6))
-    info("weighted joo in nt =")
-    dump(reshape(round(gn+un, 2), 2, 6))
-    info("lambda in nt =")
-    dump(reshape(round(la, 2), 2, 6))
-    # complementarity function
-    cn = 1.0
-#   C = la - clamp(la - cn*(gn+un), 0, Inf)
-#   C = la + clamp(la - cn*(gn+un), 0, Inf)
-    C = la + cn*(un - gn)
-    info("complementarity function =")
-    dump(reshape(round(C, 2), 2, 6))
-
-    g = zeros(length(gn))
-
-    for i=1:2:dim
-        if i == 1
-        #if C[i] > 0
-            if i == 7
-                info("skipping root dof 7")
-                continue
-            end
-            info("dof $i in active set")
-#           g[i] =  -gn[i]
-        else
-            info("dof $i not in active set")
-#           C1[i,:] = 0
-#           C2[i,:] = 0
-        end
-    end
-
-    g[1] = 2.0
-    g[3] = 2.0
-#=
-    for i=2:2:dim
-        C1[i,:] = 0
-        C2[i,:] = 0
-    end
-=#
-    d = [7, 8]
-    C2[d, :] = 0
-
-    assembly.g = sparse(g)
-    assembly.C1 = C1
-    assembly.C2 = C2
-    info("PDASS ready.")
-end
-
-function linear_system_solver_preprocess!(solver, iter, time, K, f, C1, C2, D, g, sol, la, ::Type{Val{:before_solution}})
-    # example how to use preprocessor to dump matrices before solution
-    @debug begin
-        info("stiffness matrix")
-        dump(round(full(K), 3))
-        info("constraint matrix C1")
-        dump(round(full(C1), 3))
-        info("constraint matrix C2")
-        dump(round(full(C2), 3))
-        info("force vector")
-        dump(round(full(f)', 3))
-        info("constraint vector")
-        dump(round(full(g)', 3))
-    end
-end
-
-function linear_system_solver_postprocess!(solver, iter, time, K, f, C1, C2, D, g, x, la, ::Type{Val{:after_solution}})
-    @debug begin
-        info("solution vector")
-        dump(round(full(x)', 3))
-        info("reaction force vector")
-        dump(round(full(la)', 3))
-    end
-end
-
-@testset "2d frictionless contact" begin
-    gap = [1.0, 0.0]
-    nodes = Node[
-        [6.0, 6.0],
-        [6.0, 12.0]+gap,
-        [0.0, 0.0],
-        [6.0, 0.0],
-        [6.0, 0.0]+gap,
-        [18.0, 0.0]+gap]
-    fel1 = Tri3([1, 3, 4])
-    fel2 = Tri3([2, 5, 6])
-    force = Seg2([3, 1])
-    bnd1 = Seg2([3, 4])
-    bnd2 = Seg2([5, 6])
-    sel = Seg2([1, 4])
-    mel = Seg2([2, 5])
-    update!([fel1, fel2, force, sel, mel], "geometry", nodes)
-    update!([bnd1, bnd2], "geometry", nodes)
-
-    prob = FieldProblem(PlaneStressLinearElasticityProblem, "bodies", 2)
-    push!(prob, fel1, fel2)
-    push!(prob, force)
-    update!([fel1, fel2], "youngs modulus", 90.0)
-    update!([fel1, fel2], "poissons ratio", 0.25)
-    update!([force], "displacement traction force 1", 2*6/sqrt(2))
-
-    bc = BoundaryProblem(DirichletProblem, "support", "displacement", 2)
-    push!(bc, bnd1, bnd2)
-    update!(get_elements(bc), "displacement", 0.0 => Vector{Float64}[[0.0, 0.0], [0.0, 0.0]])
-
-    cont = BoundaryProblem(MortarProblem, "contact", "displacement", 2)
-    push!(cont, sel, mel)
-    calculate_normal_tangential_coordinates!(sel, 0.0)
-    nt = sel("normal-tangential coordinates", [0.0], 0.0)
-    info("normal direction = $(nt)")
-    sel["master elements"] = [mel]
-    # remove coefficients from node 4 (dofs 7-8) because this conflicts with dirichlet bc.
-    # add_postprocessor!(cont, :remove_constraint_from_dofs, [7, 8])
-    # remove tangential direction constraints
-    # add_postprocessor!(cont, :remove_tangential_constraints)
-    # apply PDASS
-    add_postprocessor!(cont, :primal_dual_active_set_strategy)
-
-    @debug begin
-        info("fel1.fields = $(fel1.fields)")
-    end
-
-    solver = DirectSolver()
-    push!(solver, prob)
-    push!(solver, bc)
-    push!(solver, cont)
-    solver.solve_residual = false
-    set_linear_system_solver!(solver, :UMFPACK)
-    set_nonlinear_max_iterations!(solver, 5)
-    add_linear_system_solver_preprocessor!(solver, :before_solution)
-    add_linear_system_solver_postprocessor!(solver, :after_solution)
-    add_linear_system_solver_preprocessor!(solver, :dump_matrices)
-    time = 0.0
-    call(solver, time)
-
-    @debug begin
-        u = calculate_nodal_vector("displacement", 2, get_elements(prob), time)
-        info("solution vector")
-        dump(reshape(round(u, 8), 2, 6))
-#       la = calculate_nodal_vector("reaction force", 2, get_elements(prob), time)
-#       info("reaction force")
-#       dump(reshape(round(la, 8), 2, 6))
-    end
-end
-
diff --git a/test/test_mortar_2d_mesh_tie.jl b/test/test_mortar_2d_mesh_tie.jl
index e40c786..3c67c27 100644
--- a/test/test_mortar_2d_mesh_tie.jl
+++ b/test/test_mortar_2d_mesh_tie.jl
@@ -99,7 +99,8 @@ end
 
 
 #=
-
+# TODO: if one forget plane_stress solver gives singular exception and it's
+    hard to trace to the source of problem
 @testset "expect clear error when trying to solve 2d model in 3d setting" begin
     p1, p2, p3, p4 = get_test_model()
 #   p1.properties.formulation = :plane_stress
@@ -115,5 +116,164 @@ end
     info("u = $u")
     @test isapprox(u, [0.0, 0.05])
 end
-
 =#
+
+function JuliaFEM.get_mesh(::Type{Val{Symbol("1x1 block splitted to upper and lower")}})
+    meshfile = Pkg.dir("JuliaFEM") * "/test/testdata/block_2d.med"
+    mesh = aster_read_mesh(meshfile)
+end
+
+function JuliaFEM.get_model(::Type{Val{Symbol("splitted block, plane stress elasticity and mesh tie")}})
+    mesh = get_mesh("1x1 block splitted to upper and lower")
+
+    upper = Problem(Elasticity, "upper", 2)
+    upper.properties.formulation = :plane_stress
+    upper.elements = create_elements(mesh, "UPPER")
+    update!(upper.elements, "youngs modulus", 100.0)
+    update!(upper.elements, "poissons ratio", 1/3)
+
+    lower = Problem(Elasticity, "lower", 2)
+    lower.properties.formulation = :plane_stress
+    lower.elements = create_elements(mesh, "LOWER")
+    update!(lower.elements, "youngs modulus", 100.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.1)
+    update!(bc_upper.elements, "displacement 2", -0.1)
+
+    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)
+
+    bc_corner = Problem(Dirichlet, "fix model from lower left corner to prevent singularity", 2, "displacement")
+    node_ids = find_nearest_nodes(mesh, [0.0, 0.0])
+    bc_corner.elements = [Element(Poi1, node_ids)]
+    update!(bc_corner.elements, "geometry", mesh.nodes)
+    update!(bc_corner.elements, "displacement 1", 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")
+    update!(interface_slave_elements, "master elements", interface_master_elements)
+    interface.elements = [interface_master_elements; interface_slave_elements]
+
+    solver = Solver(Nonlinear)
+    push!(solver, upper, lower, bc_upper, bc_lower, interface, bc_corner)
+
+    return solver
+
+end
+
+@testset "test mesh tie with splitted block and plane stress elasticity" begin
+    solver = get_model("splitted block, plane stress elasticity and mesh tie")
+    upper, lower, bc_upper, bc_lower, interface = solver.problems
+    call(solver)
+    @test solver.properties.iteration == 2
+    slave_elements = get_slave_elements(interface)
+    node_ids, la = get_nodal_vector(slave_elements, "reaction force", 0.0)
+    for lai in la
+        @test isapprox(lai, [0.0, 10.0])
+    end
+end
+
+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)
+
+    mesh = get_mesh("curved 2d block splitted to upper and lower")
+
+    upper = Problem(Elasticity, "upper", 2)
+    upper.properties.formulation = :plane_stress
+    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.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.tolerance = tolerance
+    interface.properties.rotate_normals = rotate_normals
+    interface.properties.dual_basis = dual_basis
+
+    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=true, tolerance=10, dy=0.0, rotate_normals=true,
+        dual_basis=false)
+    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)
+    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)
+    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)
+    call(solver)
+    interface = solver["interface between upper and lower block"]
+    @test solver.properties.iteration == 2
+    @test isapprox(norm(interface.assembly.u), 0.34318800698017704)
+end
+
diff --git a/test/test_nodal_constraints.jl b/test/test_nodal_constraints.jl
new file mode 100644
index 0000000..f9c502c
--- /dev/null
+++ b/test/test_nodal_constraints.jl
@@ -0,0 +1,82 @@
+
+# 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.Test
+
+function JuliaFEM.get_model(::Type{Val{Symbol("1x1 plane stress quad4 block")}})
+    
+    X = Dict{Int, Vector{Float64}}(
+        1 => [0.0, 0.0],
+        2 => [1.0, 0.0],
+        3 => [1.0, 1.0],
+        4 => [0.0, 1.0])
+
+    body = Problem(Elasticity, "body", 2)
+    body.properties.formulation = :plane_stress
+    body.elements = [Element(Quad4, [1, 2, 3, 4])]
+    update!(body.elements, "geometry", X)
+    update!(body.elements, "youngs modulus", 288.0)
+    update!(body.elements, "poissons ratio", 1/3)
+
+    # boundary conditions
+    bc_13 = Problem(Dirichlet, "symmetry 13", 2, "displacement")
+    bc_13.properties.dual_basis = true
+    bc_13.elements = [Element(Seg2, [1, 2])]
+    update!(bc_13.elements, "geometry", X)
+    update!(bc_13.elements, "displacement 2", 0.0)
+
+    bc_23 = Problem(Dirichlet, "symmetry 23", 2, "displacement")
+    bc_23.properties.dual_basis = true
+    bc_23.elements = [Element(Seg2, [4, 1])]
+    update!(bc_23.elements, "geometry", X)
+    update!(bc_23, "displacement 1", 0.0)
+
+    solver = Solver(Nonlinear, "1x1 plane stress quad4 block")
+    push!(solver, body, bc_13, bc_23)
+
+    return solver
+end
+
+@testset "test dirichlet spc in point" begin
+    X = Dict{Int, Vector{Float64}}(
+        1 => [0.0, 0.0],
+        2 => [1.0, 0.0],
+        3 => [1.0, 1.0],
+        4 => [0.0, 1.0])
+    solver = get_model("1x1 plane stress quad4 block")
+    update!(solver["symmetry 13"], "displacement 1", 0.0)
+    update!(solver["symmetry 23"], "displacement 2", 0.0)
+    nodal_bc = Problem(Dirichlet, "dx=0.5", 2, "displacement")
+    nodal_bc.elements = [Element(Poi1, [3])]
+    update!(nodal_bc, "geometry", X)
+    update!(nodal_bc, "displacement 1", 0.5)
+    update!(nodal_bc, "displacement 2", 0.0)
+    push!(solver, nodal_bc)
+    call(solver)
+    pel = nodal_bc.elements[1]
+    la = pel("reaction force", [0.0], 0.0)
+    info("reaction force: $la")
+    info(solver["body"].assembly.u)
+    @test isapprox(pel("displacement", [], 0.0), [0.5, 0.0])
+end
+
+@testset "test nodal point force" begin
+    X = Dict{Int, Vector{Float64}}(
+        1 => [0.0, 0.0],
+        2 => [1.0, 0.0],
+        3 => [1.0, 1.0],
+        4 => [0.0, 1.0])
+    solver = get_model("1x1 plane stress quad4 block")
+    update!(solver["symmetry 13"], "displacement 1", 0.0)
+    update!(solver["symmetry 23"], "displacement 2", 0.0)
+    point_load = Element(Poi1, [3])
+    update!(point_load, "geometry", X)
+    update!(point_load, "displacement traction force 1", 72.0)
+    update!(point_load, "displacement traction force 2", 27.0)
+    push!(solver["body"], point_load)
+    call(solver)
+    @test isapprox(point_load("displacement", [], 0.0), [0.5, 0.0])
+end
+
diff --git a/test/test_preprocess_aster_reader.jl b/test/test_preprocess_aster_reader.jl
index a7f2dab..dba65e8 100644
--- a/test/test_preprocess_aster_reader.jl
+++ b/test/test_preprocess_aster_reader.jl
@@ -106,11 +106,46 @@ end
     @test mesh["connectivity"][2] == (:SE2, :GRP1, [3, 4])
 end
 
-@testset "test reading aster .med file" begin
-    fn = Pkg.dir("JuliaFEM")*"/geometry/2d_block/BLOCK_1elem.med"
+function JuliaFEM.get_mesh(::Type{Val{Symbol("block_2d_1elem_quad4")}})
+    fn = Pkg.dir("JuliaFEM") * "/test/testdata/block_2d_1elem_quad4.med"
     mesh = aster_read_mesh(fn)
-    @test haskey(mesh.element_sets, "BLOCK")
+    return mesh
+end
+
+@testset "test reading aster .med file" begin
+    mesh = get_mesh("block_2d_1elem_quad4")
+    info("nodes")
+    for (k, v) in mesh.nodes
+        info("$k => $v")
+    end
+    info("node sets")
+    for (k, v) in mesh.node_sets
+        info("$k => $v")
+    end
+    info("elements")
+    for (k, v) in mesh.elements
+        info("$k => $v, type = $(mesh.element_types[k])")
+    end
+    info("element sets")
+    for (k, v) in mesh.element_sets
+        info("$k => $v")
+    end
+    @test length(mesh.element_sets) == 5
+    @test length(mesh.node_sets) == 4
     @test length(mesh.elements) == 5
+    @test length(mesh.nodes) == 4
+    for elset in ["BLOCK", "TOP", "BOTTOM", "LEFT", "RIGHT"]
+        @test haskey(mesh.element_sets, elset)
+        @test length(mesh.element_sets[elset]) == 1
+    end
+    for nset in ["TOP_LEFT", "TOP_RIGHT", "BOTTOM_LEFT", "BOTTOM_RIGHT"]
+        @test haskey(mesh.node_sets, nset)
+        @test length(mesh.node_sets[nset]) == 1
+    end
+end
+
+@testset "test filter by element set" begin
+    mesh = get_mesh("block_2d_1elem_quad4")
     mesh2 = filter_by_element_set(mesh, "BLOCK")
     @test haskey(mesh2.element_sets, "BLOCK")
     @test length(mesh2.elements) == 1
diff --git a/test/test_von_mises_material.jl b/test/test_von_mises_material.jl
index c696883..7b2cabd 100644
--- a/test/test_von_mises_material.jl
+++ b/test/test_von_mises_material.jl
@@ -121,6 +121,7 @@ function test_von_mises_3D_basic()
 
     info("Calculation finished")
     #PyPlot.plot(ee, ss)
+    #=
     plot3D(eig_vals[:, 1], eig_vals[:, 2], eig_vals[:, 3], color="red")
     PyPlot.title("Stress path and von Mises yield surface")
     PyPlot.xlabel("Eig Stress 1")
@@ -128,6 +129,7 @@ function test_von_mises_3D_basic()
     PyPlot.zlabel("Eig Stress 3")
     PyPlot.grid()
     PyPlot.show()
+    =#
 end
 
 function test_von_mises_planestress_basic()
@@ -228,10 +230,12 @@ function test_von_mises_planestress_basic()
         push!(x_vals, s11)
         push!(y_vals, s22)
     end
+    #=
     PyPlot.plot(x_vals, y_vals)
     PyPlot.plot(ee, ss)
     PyPlot.grid()
     PyPlot.show()
+    =#
 end
 
 # test_von_mises_3D_basic()
diff --git a/test/testdata/2d_cylinder_roller_full.med b/test/testdata/2d_cylinder_roller_full.med
new file mode 100644
index 0000000..c932510
Binary files /dev/null and b/test/testdata/2d_cylinder_roller_full.med differ
diff --git a/test/testdata/block_2d.med b/test/testdata/block_2d.med
new file mode 100644
index 0000000..55e8b7a
Binary files /dev/null and b/test/testdata/block_2d.med differ
diff --git a/test/testdata/block_2d.png b/test/testdata/block_2d.png
new file mode 100644
index 0000000..c1cbc88
Binary files /dev/null and b/test/testdata/block_2d.png differ
diff --git a/test/testdata/block_2d_1elem_quad4.med b/test/testdata/block_2d_1elem_quad4.med
new file mode 100644
index 0000000..7948982
Binary files /dev/null and b/test/testdata/block_2d_1elem_quad4.med differ
diff --git a/test/testdata/block_2d_1elem_quad4.png b/test/testdata/block_2d_1elem_quad4.png
new file mode 100644
index 0000000..520c753
Binary files /dev/null and b/test/testdata/block_2d_1elem_quad4.png differ
diff --git a/test/testdata/block_2d_curved.med b/test/testdata/block_2d_curved.med
new file mode 100644
index 0000000..2c4927f
Binary files /dev/null and b/test/testdata/block_2d_curved.med differ
diff --git a/test/testdata/block_2d_curved_dense.med b/test/testdata/block_2d_curved_dense.med
new file mode 100644
index 0000000..8a560bf
Binary files /dev/null and b/test/testdata/block_2d_curved_dense.med differ
diff --git a/test/testdata/block_3d.med b/test/testdata/block_3d.med
new file mode 100644
index 0000000..63c0070
Binary files /dev/null and b/test/testdata/block_3d.med differ
diff --git a/test/testdata/hertz_2d_full.med b/test/testdata/hertz_2d_full.med
new file mode 100644
index 0000000..2733f0b
Binary files /dev/null and b/test/testdata/hertz_2d_full.med differ