diff --git a/src/JuliaFEM.jl b/src/JuliaFEM.jl
index 8ed736a..bfbfa35 100644
--- a/src/JuliaFEM.jl
+++ b/src/JuliaFEM.jl
@@ -17,8 +17,7 @@ using Logging
 Logging.configure(level=INFO)
 
 if haskey(ENV, "JULIAFEM_LOGLEVEL")
-    ENV["JULIAFEM_LOGLEVEL"] == "INFO" && Logging.configure(level=INFO)
-    ENV["JULIAFEM_LOGLEVEL"] == "DEBUG" && Logging.configure(level=DEBUG)
+    Logging.configure(level=LogLevel(ENV["JULIAFEM_LOGLEVEL"]))
 end
 
 export info, debug
@@ -141,16 +140,6 @@ export aster_create_elements, parse_aster_med_file, is_aster_mail_keyword,
        aster_read_mesh_names, aster_read_node_sets, aster_read_nodes, RMEDFile
 end
 
-function get_mesh(mesh_name::AbstractString, args...; kwargs...)
-    return get_mesh(Val{Symbol(mesh_name)}, args...; kwargs...)
-end
-
-function get_model(model_name::AbstractString, args...; kwargs...)
-    return get_model(Val{Symbol(model_name)}, args...; kwargs...)
-end
-
-export get_mesh, get_model
-
 module Postprocess
 
 include("postprocess_utils.jl")
diff --git a/src/problems_contact.jl b/src/problems_contact.jl
index 483bfeb..bac0010 100644
--- a/src/problems_contact.jl
+++ b/src/problems_contact.jl
@@ -25,6 +25,10 @@ type Contact <: BoundaryProblem
     remove_nodes :: Vector{Int}
     always_in_contact :: Bool
     update_contact_pairing :: Bool
+    iteration :: Int
+    contact_state_in_first_iteration :: Symbol
+    alpha :: Float64
+    drop_tolerance :: Float64
     store_fields :: Vector{AbstractString}
 end
 
@@ -47,6 +51,10 @@ function Contact()
         [],      # remove these nodes always from set
         false,   # mainly for debugging, do not remove inactive nodes
         true,    # update contact pairing on each loop
+        1,       # iteration counter
+        :AUTO,   # contact state in first iteration, AUTO, INACTIVE, ACTIVE
+        0.0,     # alpha basis transform parameter
+        1.0e-9,  # drop tolerance
         default_fields)
 end
 
@@ -55,24 +63,29 @@ function get_unknown_field_name(problem::Problem{Contact})
 end
 
 function get_formulation_type(problem::Problem{Contact})
+    #=
     if problem.properties.use_forwarddiff
         return :forwarddiff
     else
         return :incremental
     end
+    =#
+    return :incremental
+    #return :forwarddiff
 end
 
 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")
+        debug("assuming dimension of mesh tie surface is $dim")
+        debug("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}
     use_forwarddiff = Val{problem.properties.use_forwarddiff}
     assemble!(problem, time, dimension, finite_sliding, friction, use_forwarddiff)
+    problem.properties.iteration += 1
 end
 
 typealias ContactElements2D Union{Seg2}
diff --git a/src/problems_contact_2d.jl b/src/problems_contact_2d.jl
index 4f9f9b9..cac52d3 100644
--- a/src/problems_contact_2d.jl
+++ b/src/problems_contact_2d.jl
@@ -1,14 +1,62 @@
 # This file is a part of JuliaFEM.
 # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
 
+function create_rotation_matrix(element::Element{Seg2}, time::Float64)
+    n = element("normal", time)
+    R = [0.0 -1.0; 1.0 0.0]
+    t1 = R'*n[1]
+    t2 = R'*n[2]
+    Q1 = [n[1] t1]
+    Q2 = [n[2] t2]
+    Z = zeros(2, 2)
+    Q = [Q1 Z; Z Q2]
+    return Q
+end
+
+function create_contact_segmentation(problem::Problem{Contact}, slave_element::Element{Seg2}, master_elements::Vector, time::Float64; deformed=false)
+    
+    result = []
+
+    x1 = slave_element("geometry", time)
+
+    if deformed
+        x1 += slave_element("displacement", time)
+    end
+
+    for master_element in master_elements
+
+        x2 = master_element("geometry", time)
+
+        if deformed
+            x2 += master_element("displacement", time)
+        end
+
+        if norm(mean(x1) - x2[1]) / norm(x1[2] - x1[1]) > problem.properties.distval
+            continue
+        end
+        if norm(mean(x1) - x2[2]) / norm(x1[2] - x1[1]) > problem.properties.distval
+            continue
+        end
+        
+        # 3.1 calculate segmentation
+        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])
+        if isapprox(l, 0.0)
+            continue # no contribution in this master element
+        end
+        push!(result, (master_element, xi1, l))
+    end
+    return result
+end
+
 """
 Frictionless 2d small sliding contact without forwarddiff.
 
 true/false flags: finite_sliding, friction, use_forwarddiff
 """
-function assemble!(problem::Problem{Contact}, time::Float64,
-                   ::Type{Val{1}}, ::Type{Val{false}},
-                   ::Type{Val{false}}, ::Type{Val{false}}; debug=false)
+function assemble!(problem::Problem{Contact}, time::Float64, ::Type{Val{1}}, ::Type{Val{false}}, ::Type{Val{false}}, ::Type{Val{false}})
 
     props = problem.properties
     field_dim = get_unknown_field_dimension(problem)
@@ -33,51 +81,22 @@ function assemble!(problem::Problem{Contact}, time::Float64,
         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
+        Q2 = create_rotation_matrix(slave_element, time)
 
-        if "element area" in props.store_fields
-            element_area = 0.0
-            for ip in get_integration_points(slave_element)
-                detJ = slave_element(ip, time, Val{:detJ})
-                w = ip.weight*detJ
-                element_area += w
-            end
-            update!(slave_element, "element area", time => element_area)
+        master_elements = slave_element("master elements", time)
+        segmentation = create_contact_segmentation(problem, slave_element, master_elements, time)
+        if length(segmentation) == 0 # no overlapping in master and slave surfaces with this slave element
+            continue
         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
-
-            if norm(mean(X1) - X2[1]) / norm(X1[2] - X1[1]) > props.distval
-                continue
-            end
-
-            if norm(mean(X1) - X2[2]) / norm(X1[2] - X1[1]) > props.distval
-                continue
-            end
-            
-            # 3.1 calculate segmentation
-            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
+        Ae = eye(nsl)
+        if props.dual_basis
             De = zeros(nsl, nsl)
             Me = zeros(nsl, nsl)
-            Ae = zeros(nsl, nsl)
-            if props.dual_basis
+            for (master_element, xi1, l) in segmentation
                 for ip in get_integration_points(slave_element, 3)
                     detJ = slave_element(ip, time, Val{:detJ})
                     w = ip.weight*detJ*l
@@ -88,14 +107,21 @@ function assemble!(problem::Problem{Contact}, time::Float64,
                     Me += w*N1*N1'
                 end
                 Ae = De*inv(Me)
-            else
-                Ae = eye(nsl)
             end
+        end
+
+        # loop all segments
+        for (master_element, xi1, l) in segmentation
+
+            nm = length(master_element)
+            X2 = master_element("geometry", time)
+            u2 = master_element("displacement", time)
+            x2 = X2 + u2
             
             # 3.3. loop integration points of one integration segment and calculate
             # local mortar matrices
-            fill!(De, 0.0)
-            fill!(Me, 0.0)
+            De = zeros(nsl, nsl)
+            Me = zeros(nsl, nsl)
             Ne = zeros(nsl, 2*nsl)
             Te = zeros(nsl, 2*nsl)
             He = zeros(nsl, 2*nsl)
@@ -117,7 +143,7 @@ function assemble!(problem::Problem{Contact}, time::Float64,
                 t_s /= norm(t_s)
                 xi_m = project_from_slave_to_master(master_element, X_s, n_s, time)
                 N2 = vec(get_basis(master_element, xi_m, time))
-                X_m = N2*X2 
+                X_m = N2*X2
 
                 u_s = N1*u1
                 u_m = N2*u2
@@ -182,14 +208,34 @@ function assemble!(problem::Problem{Contact}, time::Float64,
 
     la = problem.assembly.la
     ndofs = length(la)
-#   info("contact ndofs: $ndofs")
-#   info("Rn = $Rn")
 
     C1 = sparse(problem.assembly.C1, ndofs, ndofs)
     C2 = sparse(problem.assembly.C2, ndofs, ndofs)
     D = sparse(problem.assembly.D, ndofs, ndofs)
     g = full(problem.assembly.g, ndofs, 1)
     c = full(problem.assembly.c, ndofs, 1)
+    
+    for j in S
+        dofs = [2*(j-1)+1, 2*(j-1)+2]
+        weighted_gap[j] = g[dofs]
+    end
+    
+    state = problem.properties.contact_state_in_first_iteration
+    if problem.properties.iteration == 1
+        info("First contact iteration, initial contact state = $state")
+
+        if state == :AUTO
+            avg_gap = mean([weighted_gap[j][1] for j in S])
+            std_gap = std([weighted_gap[j][1] for j in S])
+            if (avg_gap < 1.0e-12) && (std_gap < 1.0e-12)
+                state = :ACTIVE
+            else
+                state = :UNKNOWN
+            end
+            info("Average weighted gap = $avg_gap, std gap = $std_gap, automatically determined contact state = $state")
+        end
+
+    end
 
     # active / inactive node detection
     for j in S
@@ -204,7 +250,6 @@ function assemble!(problem::Problem{Contact}, time::Float64,
             contact_pressure[j] = [0.0, 0.0]
         end
 
-#       contact_pressure[j] = c[dofs]
         complementarity_condition[j] = contact_pressure[j] - weighted_gap[j]
         if complementarity_condition[j][1] < 0
             is_inactive[j] = 1
@@ -216,11 +261,24 @@ function assemble!(problem::Problem{Contact}, time::Float64,
             is_active[j] = 1
             is_slip[j] = 1
             is_stick[j] = 0
-#           _c1 = complementarity_condition[j][1]
-#           _c2 = c[dofs]
-#           _c3 = contact_pressure[j][1]
-#           _c4 = g[dofs]
-#           info("active $j: c1 = $_c1, c2 = $_c2, c3 = $_c3, c4 = $_c4")
+        end
+    end
+    
+    if (problem.properties.iteration == 1) && (state == :ACTIVE)
+        for j in S
+            is_inactive[j] = 0
+            is_active[j] = 1
+            is_slip[j] = 1
+            is_stick[j] = 0
+        end
+    end
+    
+    if (problem.properties.iteration == 1) && (state == :INACTIVE)
+        for j in S
+            is_inactive[j] = 1
+            is_active[j] = 0
+            is_slip[j] = 0
+            is_stick[j] = 0
         end
     end
 
@@ -246,34 +304,33 @@ function assemble!(problem::Problem{Contact}, time::Float64,
         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
+    debug("# | 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][1], 3)) | $(round(contact_pressure[j][1], 3)) | $(round(complementarity_condition[j][1], 3))"
+        debug(str1 * str2)
+    end
 
-    # solve variational inequality
-    
+    debug("normals: ", normals)
+
+    # solve variational inequality 
 
     # constitutive modelling in tangent direction, frictionless contact
-    #=
     for j in S
         dofs = [2*(j-1)+1, 2*(j-1)+2]
         if (is_active[j] == 1) && (is_slip[j] == 1)
-            info("$j is in active/slip, removing tangential constraint $(dofs[2])")
+            debug("$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
-    =#
 
     # remove inactive nodes from assembly
     for j in S
         dofs = [2*(j-1)+1, 2*(j-1)+2]
         if is_inactive[j] == 1
-#           info("$j is inactive, removing dofs $dofs")
+            debug("$j is inactive, removing dofs $dofs")
             C1[dofs,:] = 0.0
             C2[dofs,:] = 0.0
             D[dofs,:] = 0.0
diff --git a/src/problems_contact_2d_autodiff.jl b/src/problems_contact_2d_autodiff.jl
index 44c9c26..4de75f7 100644
--- a/src/problems_contact_2d_autodiff.jl
+++ b/src/problems_contact_2d_autodiff.jl
@@ -163,6 +163,38 @@ function assemble!(problem::Problem{Contact}, time::Float64,
             n1 = Field(Vector[normals[:,i] for i in slave_element_nodes])
             nnodes = size(slave_element, 2)
 
+            # construct dual basis
+            De = zeros(nnodes, nnodes)
+            Me = zeros(nnodes, nnodes)
+            for master_element in slave_element("master elements", time)
+
+                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
+
+                # calculate segmentation: we care only about endpoints
+                xi1a = project_from_master_to_slave(slave_element, x1, n1, x2[1])
+                xi1b = project_from_master_to_slave(slave_element, x1, n1, x2[2])
+                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
+
+                for ip in get_integration_points(slave_element, 3)
+                    # jacobian of slave element in deformed state
+                    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 = get_basis(slave_element, xi_s, time)
+                    De += w*diagm(vec(N1))
+                    Me += w*N1'*N1
+                end
+            end
+            
+            Ae = De*inv(Me)
+
             # 3. loop all master elements
             for master_element in slave_element("master elements", time)
 
@@ -183,21 +215,6 @@ function assemble!(problem::Problem{Contact}, time::Float64,
                 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, 3)
-                    # jacobian of slave element in deformed state
-                    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 = get_basis(slave_element, xi_s, time)
-                    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)
 
@@ -237,27 +254,67 @@ function assemble!(problem::Problem{Contact}, time::Float64,
         # 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")
+        state = problem.properties.contact_state_in_first_iteration
+        if problem.properties.iteration == 1
+            info("First contact iteration, initial contact state = $state")
+            if state == :AUTO
+                avg_gap = ForwardDiff.value(mean([gap[1, j] for j in S]))
+                std_gap = ForwardDiff.value(std([gap[1, j] for j in S]))
+                if (avg_gap < 1.0e-12) && (std_gap < 1.0e-12)
+                    state = :ACTIVE
+                else
+                    state = :UNKNOWN
+                end
+                info("Average weighted gap = $avg_gap, std gap = $std_gap, automatically determined contact state = $state")
+            end
+        end
 
-        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])
+        is_active = Dict{Int, Bool}()
+        condition = Dict()
 
-            if lan - gap[1, j] > 0
-#               info("set node $j active, normal direction = $(ForwardDiff.get_value(n)), tangent plane = $(ForwardDiff.get_value(t))")
+        for j in S
+            if j in props.always_in_contact
+                is_active[j] = true
+                continue
+            end
+            lan = dot(normals[:,j], la[:,j])
+            condition[j] = ForwardDiff.value(lan - gap[1, j])
+            is_active[j] = condition[j] > 0
+        end
+
+        if problem.properties.iteration == 1 && state == :ACTIVE
+            for j in S
+                is_active[j] = true
+            end
+        end
+
+        if problem.properties.iteration == 1 && state == :INACTIVE
+            for j in S
+                is_active[j] = false
+            end
+        end
+
+        if Logging._root.level == DEBUG
+            debug("Summary of nodes")
+            for j in sort(collect(keys(is_active)))
+                n = map(ForwardDiff.value, normals[:,j])
+                debug("$j, c=$(condition[j]), s=$(is_active[j]), n=$n")
+            end
+        end
+
+        for j in S
+
+            if is_active[j]
+                n = normals[:,j]
+                t = Q'*n
+                lan = dot(n, la[:,j])
+                lat = dot(t, la[:,j])
                 C[1,j] += gap[1, j]
                 C[2,j] += lat
             else
                 C[:,j] = la[:,j]
             end
+
         end
 
         return vec([fc C])
@@ -269,6 +326,7 @@ function assemble!(problem::Problem{Contact}, time::Float64,
     if length(x) == 0
         error("2d autodiff contact problem: initialize problem.assembly.u & la before solution")
     end
+
     A = ForwardDiff.jacobian(calculate_interface, x)
     b = calculate_interface(x)
     A = sparse(A)
@@ -278,19 +336,40 @@ function assemble!(problem::Problem{Contact}, time::Float64,
 
     ndofs = round(Int, length(x)/2)
     K = A[1:ndofs,1:ndofs]
-    C1 = transpose(A[1:ndofs,ndofs+1:end])
+    C1 = 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]
+    
+    f += C1*problem.assembly.la
+    g += D*problem.assembly.la
 
-    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)
+    #=
+    if !haskey(problem, "contact force")
+        problem.fields["contact force"] = Field(time => f)
+    else
+        update!(problem.fields["contact force"], time => f)
+    end
+
+    fc = problem.fields["contact force"]
+
+    if length(fc) > 1
+        # kick in generalized alpha rule for time integration
+        alpha = 0.5
+        info("Applying Generalized alpha time integration")
+        K = (1-alpha)*K
+        C1 = (1-alpha)*C1
+        f = alpha*fc[end-1].data
+    end
+    =#
+
+    problem.assembly.K = K
+    problem.assembly.C1 = transpose(C1)
+    problem.assembly.C2 = C2
+    problem.assembly.D = D
+    problem.assembly.f = sparse(f)
+    problem.assembly.g = sparse(g)
 
 end
 
diff --git a/src/problems_contact_3d.jl b/src/problems_contact_3d.jl
index 8bcdf12..191b790 100644
--- a/src/problems_contact_3d.jl
+++ b/src/problems_contact_3d.jl
@@ -11,6 +11,459 @@ function create_orthogonal_basis(n)
     return t1, t2
 end
 
+""" Create rotation matrix Q for element nodes rotating quantities to nt coordinaet system. """
+function create_rotation_matrix(element::Element{Tri3}, time::Float64)
+    n = element("normal", time)
+    t11, t21 = create_orthogonal_basis(n[1])
+    t12, t22 = create_orthogonal_basis(n[2])
+    t13, t23 = create_orthogonal_basis(n[3])
+    Q1_ = [n[1] t11 t21]
+    Q2_ = [n[2] t12 t22]
+    Q3_ = [n[3] t13 t23]
+    Z = zeros(3, 3)
+    Q = [
+        Q1_ Z Z
+        Z Q2_ Z
+        Z Z Q3_]
+    return Q
+end
+
+function create_rotation_matrix(element::Element{Quad4}, time::Float64)
+    n = element("normal", time)
+    t11, t21 = create_orthogonal_basis(n[1])
+    t12, t22 = create_orthogonal_basis(n[2])
+    t13, t23 = create_orthogonal_basis(n[3])
+    t14, t24 = create_orthogonal_basis(n[4])
+    Q1_ = [n[1] t11 t21]
+    Q2_ = [n[2] t12 t22]
+    Q3_ = [n[3] t13 t23]
+    Q4_ = [n[4] t14 t24]
+    Z = zeros(3, 3)
+    Q = [
+        Q1_ Z Z Z
+        Z Q2_ Z Z
+        Z Z Q3_ Z
+        Z Z Z Q4_]
+    return Q
+end
+
+function create_rotation_matrix(element::Element{Tri6}, time::Float64)
+    n = element("normal", time)
+    t11, t21 = create_orthogonal_basis(n[1])
+    t12, t22 = create_orthogonal_basis(n[2])
+    t13, t23 = create_orthogonal_basis(n[3])
+    t14, t24 = create_orthogonal_basis(n[4])
+    t15, t25 = create_orthogonal_basis(n[5])
+    t16, t26 = create_orthogonal_basis(n[6])
+    Q1_ = [n[1] t11 t21]
+    Q2_ = [n[2] t12 t22]
+    Q3_ = [n[3] t13 t23]
+    Q4_ = [n[4] t14 t24]
+    Q5_ = [n[5] t15 t25]
+    Q6_ = [n[6] t16 t26]
+    Z = zeros(3, 3)
+    Q = [
+        Q1_ Z Z Z Z Z
+        Z Q2_ Z Z Z Z
+        Z Z Q3_ Z Z Z
+        Z Z Z Q4_ Z Z
+        Z Z Z Z Q5_ Z
+        Z Z Z Z Z Q6_]
+    return Q
+end
+
+""" Create a contact segmentation between one slave element and list of master elements.
+
+Returns
+-------
+
+Vector with tuples: (master_element, polygon_clip_vertices, polygon_clip_centroid, polygon_clip_area)
+"""
+function create_contact_segmentation(slave_element, master_elements, x0, n0, time::Float64; deformed=false)
+    result = []
+    x1 = slave_element("geometry", time)
+    if deformed
+        x1 += slave_element("displacement", time)
+    end
+    S = Vector[project_vertex_to_auxiliary_plane(p, x0, n0) for p in x1]
+    for master_element in master_elements
+        x2 = master_element("geometry", time)
+        if deformed
+            x2 += master_element("displacement", time)
+        end
+        M = Vector[project_vertex_to_auxiliary_plane(p, x0, n0) for p in x2]
+        P = get_polygon_clip(S, M, n0)
+        length(P) < 3 && continue # no clipping or shared edge (no volume)
+        check_orientation!(P, n0)
+        N_P = length(P)
+        P_area = sum([norm(1/2*cross(P[i]-P[1], P[mod(i,N_P)+1]-P[1])) for i=2:N_P])
+        if isapprox(P_area, 0.0)
+            error("Polygon P has zero area")
+        end
+        C0 = calculate_centroid(P)
+        push!(result, (master_element, P, C0, P_area))
+    end
+    return result
+end
+
+"Assemble linear surface element to contact problem. """
+function assemble!(problem::Problem{Contact}, slave_element::Element{Tri3}, time::Float64)
+
+    props = problem.properties
+    field_dim = get_unknown_field_dimension(problem)
+
+    nsl = length(slave_element)
+    X1 = slave_element("geometry", time)
+    u1 = slave_element("displacement", time)
+    x1 = X1 + u1
+    n1 = slave_element("normal", time)
+    la = slave_element("reaction force", time)
+
+    Q3 = create_rotation_matrix(slave_element, time)
+
+    # project slave nodes to auxiliary plane (x0, Q)
+    xi = mean(get_reference_coordinates(slave_element))
+    N = vec(get_basis(slave_element, xi, time))
+    x0 = N*X1
+    n0 = N*n1
+
+    # create contact segmentation
+    segmentation = create_contact_segmentation(slave_element, slave_element("master elements", time), x0, n0, time)
+
+    if length(segmentation) == 0 # no overlapping surface in slave and maters
+        return
+    end
+
+    Ae = eye(nsl)
+
+    if problem.properties.dual_basis # construct dual basis
+
+        De = zeros(nsl, nsl)
+        Me = zeros(nsl, nsl)
+
+        # loop all polygons
+        for (master_element, P, C0, P_area) in segmentation
+
+            # loop integration cells
+            for cell in get_cells(P, C0)
+                virtual_element = Element(Tri3, Int[])
+                update!(virtual_element, "geometry", cell)
+                for ip in get_integration_points(virtual_element, 3)
+                    detJ = virtual_element(ip, time, Val{:detJ})
+                    w = ip.weight*detJ
+                    x_gauss = virtual_element("geometry", ip, time)
+                    xi_s, alpha = project_vertex_to_surface(x_gauss, x0, n0, slave_element, X1, time)
+                    N1 = slave_element(xi_s, time)
+                    De += w*diagm(vec(N1))
+                    Me += w*N1'*N1
+                end # integration points done
+
+            end # integration cells done
+
+        end # master elements done
+
+        Ae = De*inv(Me)
+    
+        debug("Dual basis coeffients = $Ae")
+    end
+
+    # loop all polygons
+    for (master_element, P, C0, P_area) in segmentation
+
+        nm = length(master_element)
+        X2 = master_element("geometry", time)
+        u2 = master_element("displacement", time)
+        x2 = X2 + u2
+
+        De = zeros(nsl, nsl)
+        Me = zeros(nsl, nm)
+        ce = zeros(field_dim*nsl)
+        ge = zeros(field_dim*nsl)
+
+        # loop integration cells
+        for cell in get_cells(P, C0)
+            virtual_element = Element(Tri3, Int[])
+            update!(virtual_element, "geometry", cell)
+            # loop integration point of integration cell
+            for ip in get_integration_points(virtual_element, 3)
+
+                # project gauss point from auxiliary plane to master and slave element
+                x_gauss = virtual_element("geometry", ip, time)
+                xi_s, alpha = project_vertex_to_surface(x_gauss, x0, n0, slave_element, X1, time)
+                xi_m, alpha = project_vertex_to_surface(x_gauss, x0, n0, master_element, X2, time)
+
+                detJ = virtual_element(ip, time, Val{:detJ})
+                w = ip.weight*detJ
+            
+                # add contributions
+                N1 = vec(get_basis(slave_element, xi_s, time))
+                N2 = vec(get_basis(master_element, xi_m, time))
+                Phi = Ae*N1
+                De += w*Phi*N1'
+                Me += w*Phi*N2'
+                
+                x_s = N1*(X1+u1)
+                x_m = N2*(X2+u2)
+                ge += w*vec((x_m-x_s)*Phi')
+            
+            end # integration points done
+
+        end # integration cells done
+
+        # add contribution to contact virtual work
+        sdofs = get_gdofs(problem, slave_element)
+        mdofs = get_gdofs(problem, master_element)
+        nsldofs = length(sdofs)
+        nmdofs = length(mdofs)
+        D3 = zeros(nsldofs, nsldofs)
+        M3 = zeros(nsldofs, nmdofs)
+        for i=1:field_dim
+            D3[i:field_dim:end, i:field_dim:end] += De
+            M3[i:field_dim:end, i:field_dim:end] += Me
+        end
+        
+        add!(problem.assembly.C1, sdofs, sdofs, D3)
+        add!(problem.assembly.C1, sdofs, mdofs, -M3)
+        add!(problem.assembly.C2, sdofs, sdofs, Q3'*D3)
+        add!(problem.assembly.C2, sdofs, mdofs, -Q3'*M3)
+        add!(problem.assembly.g, sdofs, Q3'*ge)
+
+    end # master elements done
+
+end
+
+
+""" Assemble quadratic surface element to contact problem. """
+function assemble!(problem::Problem{Contact}, slave_element::Element{Tri6}, time::Float64)
+
+    props = problem.properties
+    field_dim = get_unknown_field_dimension(problem)
+
+    alp = props.alpha
+
+    if alp != 0.0
+        T = [
+                1.0 0.0 0.0 0.0 0.0 0.0
+                0.0 1.0 0.0 0.0 0.0 0.0
+                0.0 0.0 1.0 0.0 0.0 0.0
+                alp alp 0.0 1.0-2*alp 0.0 0.0
+                0.0 alp alp 0.0 1.0-2*alp 0.0
+                alp 0.0 alp 0.0 0.0 1.0-2*alp
+            ]
+    else
+        T = eye(6)
+    end
+    
+    nsl = length(slave_element)
+    Xs = slave_element("geometry", time)
+    n1 = slave_element("normal", time)
+
+    Q3 = create_rotation_matrix(slave_element, time)
+
+    Ae = eye(nsl)
+
+    if problem.properties.dual_basis # construct dual basis
+
+        nsl = length(slave_element)
+        De = zeros(nsl, nsl)
+        Me = zeros(nsl, nsl)
+
+        for sub_slave_element in split_quadratic_element(slave_element, time)
+
+            slave_element_nodes = get_connectivity(sub_slave_element)
+            nsl = length(sub_slave_element)
+
+            X1 = sub_slave_element("geometry", time)
+            #u1 = sub_slave_element("displacement", time)
+            #x1 = X1 + u1
+            n1 = sub_slave_element("normal", time)
+            #la = sub_slave_element("reaction force", time)
+
+            # create auxiliary plane
+            xi = mean(get_reference_coordinates(sub_slave_element))
+            N = vec(get_basis(sub_slave_element, xi, time))
+            x0 = N*X1
+            n0 = N*n1
+
+            # project slave nodes to auxiliary plane
+            S = Vector[project_vertex_to_auxiliary_plane(p, x0, n0) for p in X1]
+
+            # 3. loop all master elements
+            for master_element in slave_element("master elements", time)
+
+                Xm = master_element("geometry", time)
+
+                if norm(mean(Xs) - mean(Xm)) > problem.properties.distval
+                    continue
+                end
+
+                # split master element to linear sub-elements and loop
+                for sub_master_element in split_quadratic_element(master_element, time)
+
+                    master_element_nodes = get_connectivity(sub_master_element)
+                    nm = length(sub_master_element)
+                    X2 = sub_master_element("geometry", time)
+                    #u2 = sub_master_element("displacement", time)
+                    #x2 = X2 + u2
+
+                    # 3.1 project master nodes to auxiliary plane and create polygon clipping
+                    M = Vector[project_vertex_to_auxiliary_plane(p, x0, n0) for p in X2]
+                    P = get_polygon_clip(S, M, n0)
+                    length(P) < 3 && continue # no clipping or shared edge (no volume)
+                    check_orientation!(P, n0)
+
+                    N_P = length(P)
+                    P_area = sum([norm(1/2*cross(P[i]-P[1], P[mod(i,N_P)+1]-P[1])) for i=2:N_P])
+                    if isapprox(P_area, 0.0)
+                        error("Polygon P has zero area")
+                    end
+
+                    C0 = calculate_centroid(P)
+
+                    # 4. loop integration cells
+                    for cell in get_cells(P, C0)
+                        virtual_element = Element(Tri3, Int[])
+                        update!(virtual_element, "geometry", cell)
+                        for ip in get_integration_points(virtual_element, 3)
+                            detJ = virtual_element(ip, time, Val{:detJ})
+                            w = ip.weight*detJ
+                            x_gauss = virtual_element("geometry", ip, time)
+                            xi_s, alpha = project_vertex_to_surface(x_gauss, x0, n0, slave_element, Xs, time)
+                            N1 = vec(slave_element(xi_s, time)*T)
+                            De += w*diagm(N1)
+                            Me += w*N1*N1'
+                        end # integration points done
+
+                    end # integration cells done
+
+                end # sub master elements done
+
+            end # master elements done
+
+        end # sub slave elements done
+
+        Ae = De*inv(Me)
+    
+        debug("Dual basis coeffients = $Ae")
+    end
+
+    # split slave element to linear sub-elements and loop
+    for sub_slave_element in split_quadratic_element(slave_element, time)
+
+        slave_element_nodes = get_connectivity(sub_slave_element)
+        nsl = length(sub_slave_element)
+        X1 = sub_slave_element("geometry", time)
+        n1 = sub_slave_element("normal", time)
+            
+        # create auxiliary plane
+        xi = mean(get_reference_coordinates(sub_slave_element))
+        N = vec(get_basis(sub_slave_element, xi, time))
+        x0 = N*X1
+        n0 = N*n1
+            
+        # project slave nodes to auxiliary plane
+        S = Vector[project_vertex_to_auxiliary_plane(p, x0, n0) for p in X1]
+
+        # 3. loop all master elements
+        for master_element in slave_element("master elements", time)
+
+            Xm = master_element("geometry", time)
+
+            if norm(mean(Xs) - mean(Xm)) > problem.properties.distval
+                continue
+            end
+
+            # split master element to linear sub-elements and loop
+            for sub_master_element in split_quadratic_element(master_element, time)
+
+                master_element_nodes = get_connectivity(sub_master_element)
+                nm = length(master_element)
+                X2 = sub_master_element("geometry", time)
+                #u2 = master_element("displacement", time)
+                #x2 = X2 + u2
+
+                # 3.1 project master nodes to auxiliary plane and create polygon clipping
+                M = Vector[project_vertex_to_auxiliary_plane(p, x0, n0) for p in X2]
+                P = get_polygon_clip(S, M, n0)
+                length(P) < 3 && continue # no clipping or shared edge (no volume)
+                check_orientation!(P, n0)
+
+                N_P = length(P)
+                P_area = sum([norm(1/2*cross(P[i]-P[1], P[mod(i,N_P)+1]-P[1])) for i=2:N_P])
+                if isapprox(P_area, 0.0)
+                    error("Polygon P has zero area")
+                end
+
+                C0 = calculate_centroid(P)
+
+                # integration is done in quadratic elements
+                nsl = length(slave_element)
+                nm = length(master_element)
+                De = zeros(nsl, nsl)
+                Me = zeros(nsl, nm)
+                ge = zeros(field_dim*nsl)
+
+                # 4. loop integration cells
+                for cell in get_cells(P, C0)
+                    virtual_element = Element(Tri3, Int[])
+                    update!(virtual_element, "geometry", cell)
+
+                    # 5. loop integration point of integration cell
+                    for ip in get_integration_points(virtual_element, 3)
+
+                        # project gauss point from auxiliary plane to master and slave element
+                        x_gauss = virtual_element("geometry", ip, time)
+                        xi_s, alpha = project_vertex_to_surface(x_gauss, x0, n0, slave_element, Xs, time)
+                        xi_m, alpha = project_vertex_to_surface(x_gauss, x0, n0, master_element, Xm, time)
+
+                        detJ = virtual_element(ip, time, Val{:detJ})
+                        w = ip.weight*detJ
+            
+                        # add contributions
+                        N1 = vec(get_basis(slave_element, xi_s, time)*T)
+                        N2 = vec(get_basis(master_element, xi_m, time))
+                        Phi = Ae*N1
+
+                        De += w*Phi*N1'
+                        Me += w*Phi*N2'
+               
+                        us = slave_element("displacement", time)
+                        um = master_element("displacement", time)
+                        xs = N1*(Xs+us)
+                        xm = N2*(Xs+um)
+                        ge += w*vec((xm-xs)*Phi')
+            
+                    end # integration points done
+
+                end # integration cells done
+
+                # 6. add contribution to contact virtual work
+                sdofs = get_gdofs(problem, slave_element)
+                mdofs = get_gdofs(problem, master_element)
+                nsldofs = length(sdofs)
+                nmdofs = length(mdofs)
+                D3 = zeros(nsldofs, nsldofs)
+                M3 = zeros(nsldofs, nmdofs)
+                for i=1:field_dim
+                    D3[i:field_dim:end, i:field_dim:end] += De
+                    M3[i:field_dim:end, i:field_dim:end] += Me
+                end
+                
+                add!(problem.assembly.C1, sdofs, sdofs, D3)
+                add!(problem.assembly.C1, sdofs, mdofs, -M3)
+                add!(problem.assembly.C2, sdofs, sdofs, Q3'*D3)
+                add!(problem.assembly.C2, sdofs, mdofs, -Q3'*M3)
+                add!(problem.assembly.g, sdofs, Q3'*ge)
+
+            end # sub master elements done
+
+        end # master elements done
+
+    end # sub slave elements done
+
+end
+
+
 """
 Frictionless 3d small sliding contact.
 
@@ -21,9 +474,7 @@ finite_sliding
 friction
 use_forwarddiff
 """
-function assemble!(problem::Problem{Contact}, time::Float64,
-                   ::Type{Val{2}}, ::Type{Val{false}},
-                   ::Type{Val{false}}, ::Type{Val{false}}; debug=true)
+function assemble!(problem::Problem{Contact}, time::Float64, ::Type{Val{2}}, ::Type{Val{false}}, ::Type{Val{false}}, ::Type{Val{false}})
 
     props = problem.properties
     field_dim = get_unknown_field_dimension(problem)
@@ -36,153 +487,8 @@ function assemble!(problem::Problem{Contact}, time::Float64,
     update!(slave_elements, "normal", time => normals)
 
     # 2. loop all slave elements
-    for (slave_num, slave_element) in enumerate(slave_elements)
-
-        nsl = length(slave_element)
-        X1 = slave_element("geometry", time)
-        u1 = slave_element("displacement", time)
-        la = slave_element("reaction force", time)
-        n1 = slave_element("normal", time)
-        if nsl == 3
-            t11, t21 = create_orthogonal_basis(n1[1])
-            t12, t22 = create_orthogonal_basis(n1[2])
-            t13, t23 = create_orthogonal_basis(n1[3])
-            Q1_ = [n1[1] t11 t21]
-            Q2_ = [n1[2] t12 t22]
-            Q3_ = [n1[3] t13 t23]
-            Z = zeros(3, 3)
-            Q3 = [Q1_ Z Z; Z Q2_ Z; Z Z Q3_]
-        elseif nsl == 4
-            t11, t21 = create_orthogonal_basis(n1[1])
-            t12, t22 = create_orthogonal_basis(n1[2])
-            t13, t23 = create_orthogonal_basis(n1[3])
-            t14, t24 = create_orthogonal_basis(n1[4])
-            Q1_ = [n1[1] t11 t21]
-            Q2_ = [n1[2] t12 t22]
-            Q3_ = [n1[3] t13 t23]
-            Q4_ = [n1[4] t14 t24]
-            Z = zeros(3, 3)
-            Q3 = [Q1_ Z Z Z; Z Q2_ Z Z; Z Z Q3_ Z; Z Z Z Q4_]
-        else
-            error("nsl = $nsl")
-        end
-        contact_area = 0.0
-        contact_error = 0.0
-
-        element_area = 0.0
-        for ip in get_integration_points(slave_element)
-            detJ = slave_element(ip, time, Val{:detJ})
-            w = ip.weight*detJ
-            element_area += w
-        end
-        if "element area" in props.store_fields
-            update!(slave_element, "element area", time => element_area)
-        end
-
-#       if slave_num == 1
-#           info("First slave element area = $element_area")
-#           info("NT basis of first slave element")
-#           dump(Q3)
-#       end
-
-        # project slave nodes to auxiliary plane (x0, Q)
-        xi = mean(get_reference_coordinates(slave_element))
-        N = vec(get_basis(slave_element, xi, time))
-        x0 = N*X1
-        n0 = N*n1
-        S = Vector[project_vertex_to_auxiliary_plane(p, x0, n0) for p in X1]
-
-        # 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
-            norm(mean(X1) - X2[3]) / norm(X1[2] - X1[1]) < props.distval || continue
-            =#
-
-            # 3.1 project master nodes to auxiliary plane and create polygon clipping
-            M = Vector[project_vertex_to_auxiliary_plane(p, x0, n0) for p in X2]
-            P = get_polygon_clip(S, M, n0)
-            length(P) < 3 && continue # no clipping or shared edge (no volume)
-            check_orientation!(P, n0)
-            C0 = calculate_centroid(P)
-
-            De = zeros(nsl, nsl)
-            Me = zeros(nsl, nm)
-            ge = zeros(field_dim*nsl)
-
-            # 4. loop integration cells
-            for cell in get_cells(P, C0)
-                virtual_element = Element(Tri3, Int[])
-                update!(virtual_element, "geometry", cell)
-
-                # 5. loop integration point of integration cell
-                for ip in get_integration_points(virtual_element, 3)
-
-                    # project gauss point from auxiliary plane to master and slave element
-                    x_gauss = virtual_element("geometry", ip, time)
-                    if isnan(x_gauss[1])
-                        info("is nan")
-                        info("x_gauss = $x_gauss")
-                        info("cell = $cell")
-                        info("C0 = $C0")
-                        info("P = $P")
-                        info("S = $S")
-                        info("M = $M")
-                        info("n0 = $n0")
-                        error("nan, unable to continue")
-                    end
-                    xi_s, alpha = project_vertex_to_surface(x_gauss, x0, n0, slave_element, X1, time)
-                    xi_m, alpha = project_vertex_to_surface(x_gauss, x0, n0, master_element, X2, time)
-
-                    detJ = virtual_element(ip, time, Val{:detJ})
-                    w = ip.weight*detJ
-                    # add contributions
-                    N1 = vec(get_basis(slave_element, xi_s, time))
-                    N2 = vec(get_basis(master_element, xi_m, time))
-                    De += w*N1*N1'
-                    Me += w*N1*N2'
-                    
-                    x_s = N1*(X1+u1)
-                    x_m = N2*(X2+u2)
-                    ge += w*vec((x_m-x_s)*N1')
-                    contact_area += w
-                    n_s = N1*n1
-                    contact_error += 1/2*w*dot(n_s, x_s-x_m)^2
-                end # integration points done
-
-            end # integration cells done
-
-            # 6. add contribution to contact virtual work
-            sdofs = get_gdofs(problem, slave_element)
-            mdofs = get_gdofs(problem, master_element)
-            nsldofs = length(sdofs)
-            nmdofs = length(mdofs)
-            D3 = zeros(nsldofs, nsldofs)
-            M3 = zeros(nsldofs, nmdofs)
-            for i=1:field_dim
-                D3[i:field_dim:end, i:field_dim:end] += De
-                M3[i:field_dim:end, i:field_dim:end] += Me
-            end
-            
-            add!(problem.assembly.C1, sdofs, sdofs, D3)
-            add!(problem.assembly.C1, sdofs, mdofs, -M3)
-            add!(problem.assembly.C2, sdofs, sdofs, Q3'*D3)
-            add!(problem.assembly.C2, sdofs, mdofs, -Q3'*M3)
-            add!(problem.assembly.g, sdofs, Q3'*ge)
-
-        end # master elements done
-
-        if "contact area" in props.store_fields
-            update!(slave_element, "contact area", time => contact_area)
-        end
-
+    for slave_element in slave_elements
+        assemble!(problem, slave_element, time)
     end # slave elements done, contact virtual work ready
 
     S = sort(collect(keys(normals))) # slave element nodes
@@ -196,13 +502,87 @@ function assemble!(problem::Problem{Contact}, time::Float64,
 
     la = problem.assembly.la
     ndofs = length(la)
-    
+
     C1 = sparse(problem.assembly.C1, ndofs, ndofs)
     C2 = sparse(problem.assembly.C2, ndofs, ndofs)
     D = sparse(problem.assembly.D, ndofs, ndofs)
     g = full(problem.assembly.g, ndofs, 1)
     c = full(problem.assembly.c, ndofs, 1)
 
+    maxdim = maximum(size(C1))
+    if problem.properties.alpha != 0.0
+        debug("mortar_3d: size C1 = ", size(C1), " max dim = $maxdim")
+        debug("alpha != 0.0, applying transformation D = Dh*T^-1")
+        alp = problem.properties.alpha
+        Te = [
+                1.0 0.0 0.0 0.0 0.0 0.0
+                0.0 1.0 0.0 0.0 0.0 0.0
+                0.0 0.0 1.0 0.0 0.0 0.0
+                alp alp 0.0 1.0-2*alp 0.0 0.0
+                0.0 alp alp 0.0 1.0-2*alp 0.0
+                alp 0.0 alp 0.0 0.0 1.0-2*alp
+            ]
+        invTe = [
+            1.0 0.0 0.0 0.0 0.0 0.0
+            0.0 1.0 0.0 0.0 0.0 0.0
+            0.0 0.0 1.0 0.0 0.0 0.0
+            -alp/(1-2*alp) -alp/(1-2*alp) 0.0 1/(1-2*alp) 0.0 0.0
+            0.0 -alp/(1-2*alp) -alp/(1-2*alp) 0.0 1/(1-2*alp) 0.0
+            -alp/(1-2*alp) 0.0 -alp/(1-2*alp) 0.0 0.0 1/(1-2*alp)
+        ]
+        # construct global transformation matrices T and invT
+        T = SparseMatrixCOO()
+        invT = SparseMatrixCOO()
+        for element in slave_elements
+            dofs = get_gdofs(problem, element)
+            for i=1:field_dim
+                ldofs = dofs[i:field_dim:end]
+                add!(T, ldofs, ldofs, Te)
+                add!(invT, ldofs, ldofs, invTe)
+            end
+        end
+        T = sparse(T, maxdim, maxdim, (a, b) -> b)
+        invT = sparse(invT, maxdim, maxdim, (a, b) -> b)
+        # fill diagonal
+        d = ones(size(T, 1))
+        d[get_nonzero_rows(T)] = 0.0
+        T += spdiagm(d)
+        invT += spdiagm(d)
+        #invT2 = sparse(inv(full(T)))
+        #info("invT == invT2? ", invT == invT2)
+        #maxabsdiff = maximum(abs(invT - invT2))
+        #info("max diff = $maxabsdiff")
+        C1 = C1*invT
+        C2 = C2*invT
+    end
+
+    tol = problem.properties.drop_tolerance
+    debug("Dropping small values from C1 & C2, tolerace = $tol")
+    SparseArrays.droptol!(C1, tol)
+    SparseArrays.droptol!(C2, tol)
+
+    for j in S
+        dofs = [3*(j-1)+1, 3*(j-1)+2, 3*(j-1)+3]
+        weighted_gap[j] = g[dofs]
+    end
+
+    state = problem.properties.contact_state_in_first_iteration
+    if problem.properties.iteration == 1
+        info("First contact iteration, initial contact state = $state")
+
+        if state == :AUTO
+            avg_gap = mean([weighted_gap[j][1] for j in S])
+            std_gap = std([weighted_gap[j][1] for j in S])
+            if (avg_gap < 1.0e-12) && (std_gap < 1.0e-12)
+                state = :ACTIVE
+            else
+                state = :UNKNOWN
+            end
+            info("Average weighted gap = $avg_gap, std gap = $std_gap, automatically determined contact state = $state")
+        end
+
+    end
+
     # active / inactive node detection
     for j in S
         dofs = [3*(j-1)+1, 3*(j-1)+2, 3*(j-1)+3]
@@ -218,7 +598,8 @@ function assemble!(problem::Problem{Contact}, time::Float64,
             contact_pressure[j] = [0.0, 0.0, 0.0]
         end
         complementarity_condition[j] = contact_pressure[j] - weighted_gap[j]
-        if complementarity_condition[j][1] < 0
+        
+        if complementarity_condition[j][1] < 0.0
             is_inactive[j] = 1
             is_active[j] = 0
             is_slip[j] = 0
@@ -231,45 +612,49 @@ function assemble!(problem::Problem{Contact}, time::Float64,
         end
     end
 
-    if "weighted gap" in props.store_fields
-        update!(slave_elements, "weighted gap", time => weighted_gap)
+    if (problem.properties.iteration == 1) && (state == :ACTIVE)
+        for j in S
+            is_inactive[j] = 0
+            is_active[j] = 1
+            is_slip[j] = 1
+            is_stick[j] = 0
+        end
     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)
+    
+    if (problem.properties.iteration == 1) && (state == :INACTIVE)
+        for j in S
+            is_inactive[j] = 1
+            is_active[j] = 0
+            is_slip[j] = 0
+            is_stick[j] = 0
+        end
     end
 
-    #=
-    info("# | active | inactive | stick | slip | gap | pres | comp")
+    info("# | active | 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))"
+        str1 = "$j | $(is_active[j]) | $(is_stick[j]) | $(is_slip[j]) | "
+        str2 = "$(round(weighted_gap[j][1], 3)) | $(round(contact_pressure[j][1], 3)) | $(round(complementarity_condition[j][1], 3))"
         info(str1 * str2)
     end
-    =#
-
-    # solve variational inequality
-
+    
+    # remove inactive nodes from assembly
+    for j in S
+        dofs = [3*(j-1)+1, 3*(j-1)+2, 3*(j-1)+3]
+        if is_inactive[j] == 1
+            debug("$j is inactive, removing dofs $dofs")
+            C1[dofs,:] = 0.0
+            C2[dofs,:] = 0.0
+            D[dofs,:] = 0.0
+            g[dofs,:] = 0.0
+        end
+    end
+    
     # constitutive modelling in tangent direction, frictionless contact
     for j in S
         dofs = [3*(j-1)+1, 3*(j-1)+2, 3*(j-1)+3]
         tdofs = dofs[[2,3]]
         if (is_active[j] == 1) && (is_slip[j] == 1)
-#           info("$j is in active/slip, removing tangential constraints $tdofs")
+            debug("$j is in active/slip, removing tangential constraints $tdofs")
             C2[tdofs,:] = 0.0
             g[tdofs] = 0.0
             normal = normals[j]
@@ -279,22 +664,9 @@ function assemble!(problem::Problem{Contact}, time::Float64,
         end
     end
 
-    # remove inactive nodes from assembly
-    for j in S
-        dofs = [3*(j-1)+1, 3*(j-1)+2, 3*(j-1)+3]
-        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
-        end
-    end
-
     problem.assembly.C1 = C1
     problem.assembly.C2 = C2
     problem.assembly.D = D
     problem.assembly.g = g
 
 end
-
diff --git a/src/solvers.jl b/src/solvers.jl
index 5781eba..94edf53 100644
--- a/src/solvers.jl
+++ b/src/solvers.jl
@@ -13,13 +13,15 @@ type Solver{S<:AbstractSolver}
     initialized :: Bool
     u :: Vector{Float64}
     la :: Vector{Float64}
+    alpha :: Float64             # generalized alpha time integration coefficient
+    fields :: Dict{AbstractString, Field}
     properties :: S
 end
 
 
 function Solver{S<:AbstractSolver}(::Type{S}, name="solver", properties...)
     variant = S(properties...)
-    solver = Solver{S}(name, 0.0, [], [], 0, nothing, false, [], [], variant)
+    solver = Solver{S}(name, 0.0, [], [], 0, nothing, false, [], [], 0.0, Dict(), variant)
     return solver
 end
 
@@ -33,11 +35,11 @@ function get_problems(solver::Solver)
     return solver.problems
 end
 
-function push!(solver::Solver, problem)
+function push!(solver::Solver, problem::Problem)
     push!(solver.problems, problem)
 end
 
-function getindex(solver::Solver, problem_name)
+function getindex(solver::Solver, problem_name::String)
     for problem in get_problems(solver)
         if problem.name == problem_name
             return problem
@@ -46,6 +48,10 @@ function getindex(solver::Solver, problem_name)
     throw(KeyError(problem_name))
 end
 
+function haskey(solver::Solver, field_name::String)
+    return haskey(solver.fields, field_name)
+end
+
 # one-liner helpers to identify problem types
 
 is_field_problem(problem) = false
@@ -312,6 +318,23 @@ function solve!(solver::Solver; empty_assemblies_before_solution=true, symmetric
         end
         gc()
     end
+    
+    if !haskey(solver, "fint")
+        solver.fields["fint"] = Field(time => f)
+    else
+        update!(solver.fields["fint"], time => f)
+    end
+
+    fint = solver.fields["fint"]
+
+    if length(fint) > 1
+        # kick in generalized alpha rule for time integration
+        alpha = solver.alpha
+        debug("Using generalized-α time integration, α=$alpha")
+        K = (1-alpha)*K
+        C1 = (1-alpha)*C1
+        f = (1-alpha)*f + alpha*fint[end-1].data
+    end
 
     ndofs = solver.ndofs
     u = zeros(ndofs)
diff --git a/test/test_contact_2d_small_sliding.jl b/test/test_contact_2d_small_sliding.jl
deleted file mode 100644
index 484bd8e..0000000
--- a/test/test_contact_2d_small_sliding.jl
+++ /dev/null
@@ -1,148 +0,0 @@
-# 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.Testing
-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, "youngs modulus", 96.0)
-    update!(upper, "poissons ratio", 1/3)
-
-    lower = Problem(Elasticity, "lower", 2)
-    lower.properties.formulation = :plane_stress
-    lower.elements = create_elements(mesh, "LOWER")
-    update!(lower, "youngs modulus", 96.0)
-    update!(lower, "poissons ratio", 1/3)
-
-    bc_upper = Problem(Dirichlet, "upper boundary", 2, "displacement")
-    bc_upper.elements = create_elements(mesh, "UPPER_TOP")
-    update!(bc_upper, "displacement 1", 0.0)
-    update!(bc_upper, "displacement 2", -0.15)
-
-    bc_lower = Problem(Dirichlet, "lower boundary", 2, "displacement")
-    bc_lower.elements = create_elements(mesh, "LOWER_BOTTOM")
-    update!(bc_lower, "displacement 1", 0.0)
-    update!(bc_lower, "displacement 2", 0.0)
-
-    contact = Problem(Contact, "contact between upper and lower block", 2, "displacement")
-    contact.properties.dimension = 1
-    contact.properties.rotate_normals = true
-    contact_slave_elements = create_elements(mesh, "LOWER_TOP")
-    contact_master_elements = create_elements(mesh, "UPPER_BOTTOM")
-    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_upper, bc_lower, contact)
-    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")
-    solver()
-    upper, lower, bc_upper, bc_lower, contact = solver.problems
-    @test isapprox(norm(contact.assembly.u), 0.49563347601324315)
-end
-
-function get_model(::Type{Val{Symbol("hertz contact, full 2d model")}})
-    meshfile = Pkg.dir("JuliaFEM") * "/test/testdata/hertz_2d_full.med"
-    mesh = aster_read_mesh(meshfile)
-
-    upper = Problem(Elasticity, "CYLINDER", 2)
-    upper.properties.formulation = :plane_strain
-    upper.elements = create_elements(mesh, "CYLINDER")
-    update!(upper, "youngs modulus", 70.0e3)
-    update!(upper, "poissons ratio", 0.3)
-
-    lower = Problem(Elasticity, "BLOCK", 2)
-    lower.properties.formulation = :plane_strain
-    lower.elements = create_elements(mesh, "BLOCK")
-    update!(lower, "youngs modulus", 210.0e3)
-    update!(lower, "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, "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, "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, "displacement traction force 2", -35.0e3)
-
-    contact = Problem(Contact, "contact between block and cylinder", 2, "displacement")
-    contact.properties.rotate_normals = true
-    contact.properties.finite_sliding = false
-    contact.properties.friction = false
-    contact.properties.use_forwarddiff = false
-    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
-    # from fenet d3613 advanced finite element contact benchmarks
-    # a = 6.21 mm, pmax = 3585 MPa
-    # this is a very sparse 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
-    solver = get_model("hertz contact, full 2d model")
-    upper, lower, bc_fixed, bc_sym_23, load, contact = solver.problems
-    solver()
-    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)
-    # 12 % error in maximum pressure
-    @test isapprox(maximum(pres), 3585.0; rtol = 12.0e-2)
-    # 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
-    # under 0.15 % error in reaction force
-    @test isapprox(Rn, 35.0e3; rtol=0.15e-2)
-    @test isapprox(Rt, 0.0; atol=10.0)
-end
-
diff --git a/test/test_contact_3d_small_sliding.jl b/test/test_contact_3d_small_sliding.jl
deleted file mode 100644
index 9a8ec08..0000000
--- a/test/test_contact_3d_small_sliding.jl
+++ /dev/null
@@ -1,60 +0,0 @@
-# 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.Testing
-
-@testset "3d upper side curved contact" begin
-
-    # TODO: accurate solution is not known, verify using another fem software
-    # however results look very meaningful and probably this is right.
-
-    meshfile = Pkg.dir("JuliaFEM") * "/test/testdata/block_3d_curved.med"
-    mesh = aster_read_mesh(meshfile)
-
-    upper = Problem(Elasticity, "upper", 3)
-    upper.elements = create_elements(mesh, "UPPER")
-    update!(upper, "youngs modulus", 96.0)
-    update!(upper, "poissons ratio", 1/3)
-
-    lower = Problem(Elasticity, "lower", 3)
-    lower.elements = create_elements(mesh, "LOWER")
-    update!(lower, "youngs modulus", 96.0)
-    update!(lower, "poissons ratio", 1/3)
-
-    bc_upper = Problem(Dirichlet, "upper boundary", 3, "displacement")
-    bc_upper.elements = create_elements(mesh, "UPPER_TOP")
-    update!(bc_upper, "displacement 1",  0.0)
-    update!(bc_upper, "displacement 2",  0.0)
-    update!(bc_upper, "displacement 3", -0.1)
-
-    bc_lower = Problem(Dirichlet, "lower boundary", 3, "displacement")
-    bc_lower.elements = create_elements(mesh, "LOWER_BOTTOM")
-    update!(bc_lower, "displacement 1", 0.0)
-    update!(bc_lower, "displacement 2", 0.0)
-    update!(bc_lower, "displacement 3", 0.0)
-
-    contact = Problem(Contact, "contact between upper and lower block", 3, "displacement")
-    contact_slave_elements = create_elements(mesh, "LOWER_TOP")
-    contact_master_elements = create_elements(mesh, "UPPER_BOTTOM")
-    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_upper, bc_lower, contact)
-
-    solver()
-    for element in get_slave_elements(contact)
-        normal = element("normal", [1/3, 1/3], solver.time)
-        @test isapprox(normal, [0.0, 0.0, 1.0])
-        pres = dot(normal, element("reaction force", [1/3, 1/3], solver.time))
-        info("pressure = $pres")
-        #info(element("displacement", [1/3, 1/3], solver.time))
-    end
-    normu = norm(contact.assembly.u)
-    info("displacement field norm = $normu")
-    @test isapprox(normu, 0.7417568870648232)
-
-end
diff --git a/test/test_mortar_2d_contact.jl b/test/test_mortar_2d_contact.jl
index b58867b..83f1459 100644
--- a/test/test_mortar_2d_contact.jl
+++ b/test/test_mortar_2d_contact.jl
@@ -5,7 +5,8 @@ using JuliaFEM
 using JuliaFEM.Preprocess
 using JuliaFEM.Testing
 
-function JuliaFEM.get_mesh(::Type{Val{Symbol("two elements 1.0x0.5 with 0.1 gap in y direction")}})
+function get_model()
+    
     mesh = Mesh()
     add_node!(mesh, 1, [0.0, 0.0])
     add_node!(mesh, 2, [1.0, 0.0])
@@ -27,12 +28,6 @@ function JuliaFEM.get_mesh(::Type{Val{Symbol("two elements 1.0x0.5 with 0.1 gap
     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 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
@@ -72,7 +67,7 @@ function JuliaFEM.get_model(::Type{Val{Symbol("two element contact")}})
 end
 
 @testset "test simple two element contact" begin
-    solver = get_model("two element contact")
+    solver = get_model()
     solver()
     contact = solver["LOWER_TO_UPPER"]
     master = first(contact.elements)
@@ -82,5 +77,4 @@ end
     info("u = $u, la = $la")
     @test isapprox(u, [-0.2, -0.15])
     @test isapprox(la, [0.0, -30.375])
-    # FIXME
 end
diff --git a/test/test_mortar_2d_mesh_tie.jl b/test/test_mortar_2d_mesh_tie.jl
index 64f7d21..66ee551 100644
--- a/test/test_mortar_2d_mesh_tie.jl
+++ b/test/test_mortar_2d_mesh_tie.jl
@@ -131,13 +131,9 @@ end
 end
 =#
 
-function JuliaFEM.get_mesh(::Type{Val{Symbol("1x1 block splitted to upper and lower")}})
+@testset "test mesh tie with splitted block and plane stress elasticity" begin
     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
@@ -175,113 +171,11 @@ function JuliaFEM.get_model(::Type{Val{Symbol("splitted block, plane stress elas
 
     solver = Solver(Linear)
     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
+    
     solver()
     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
-    # FIXME
-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, use_forwarddiff=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.distval = tolerance
-    interface.properties.rotate_normals = rotate_normals
-    interface.properties.dual_basis = dual_basis
-    interface.properties.use_forwarddiff = use_forwarddiff
-
-    solver = Solver(Linear)
-    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)
-    solver()
-    interface = solver["interface between upper and lower block"]
-    @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)
-    solver()
-    interface = solver["interface between upper and lower block"]
-    @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)
-    solver()
-    interface = solver["interface between upper and lower block"]
-    @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)
-    solver()
-    interface = solver["interface between upper and lower block"]
-    @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
index 17b75b1..49d48cf 100644
--- a/test/test_mortar_2d_mesh_tie_forwarddiff.jl
+++ b/test/test_mortar_2d_mesh_tie_forwarddiff.jl
@@ -6,17 +6,13 @@ using JuliaFEM.Preprocess
 using JuliaFEM.Postprocess
 using JuliaFEM.Testing
 
-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")}};
+function 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")
+    meshfile = Pkg.dir("JuliaFEM") * "/test/testdata/block_2d_curved.med"
+    mesh = aster_read_mesh(meshfile)
 
     upper = Problem(Elasticity, "upper", 2)
     upper.properties.formulation = :plane_stress
diff --git a/test/test_preprocess.jl b/test/test_preprocess.jl
index 2509b88..a80ef49 100644
--- a/test/test_preprocess.jl
+++ b/test/test_preprocess.jl
@@ -44,3 +44,94 @@ end
     nid = find_nearest_node(mesh, [0.0, 0.5]; node_set="UPPER_BOTTOM")
     @test first(nid) == 13
 end
+
+@testset "code aster / parse nodes" begin
+    section = """
+    N9  2.0 3.0 4.0
+    COOR_3D
+    N1          0.0 0.0 0.0
+    N2          1.0 0.0 0.0
+    N3          1.0 1.0 0.0
+    N4          0.0 1.0 0.0
+    N5          0.0 0.0 1.0
+    N6          1.0 0.0 1.0
+    N7          1.0 1.0 1.0
+    N8          0.0 1.0 1.0
+    FINSF
+    absdflasdf
+    N12 3.0 4.0 5.0 6.0
+    N13 3.0 4.0 5.0
+    """
+    nodes = aster_parse_nodes(section)
+    @test nodes[1] == Float64[0.0, 0.0, 0.0]
+    @test nodes[8] == Float64[0.0, 1.0, 1.0]
+    @test length(nodes) == 8
+end
+
+@testset "test reading aster .med file" begin
+    meshfile = joinpath(datadir, "block_2d_1elem_quad4.med")
+    mesh = aster_read_mesh(meshfile)
+    @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 = aster_read_mesh(joinpath(datadir, "block_2d_1elem_quad4.med"))
+    mesh2 = filter_by_element_set(mesh, :BLOCK)
+    @test haskey(mesh2.element_sets, :BLOCK)
+    @test length(mesh2.elements) == 1
+end
+
+function calculate_volume(mesh_name, eltype)
+    mesh_file = joinpath(datadir, "primitives.med")
+    mesh = aster_read_mesh(mesh_file, mesh_name)
+    elements = create_elements(mesh; element_type=eltype)
+    V = 0.0
+    time = 0.0
+    for element in elements
+        for ip in get_integration_points(element)
+            detJ = element(ip, time, Val{:detJ})
+            detJ > 0 || warn("negative determinant for element $eltype !")
+            V += ip.weight*detJ
+        end
+    end
+    info("volume of $eltype is $V")
+    return V
+end
+
+@testset "calculate volume for 1 element models" begin
+    @test isapprox(calculate_volume("TRIANGLE_TRI3_1", :Tri3), 1/2)
+    @test isapprox(calculate_volume("TRIANGLE_TRI6_1", :Tri6), 1/2)
+    @test isapprox(calculate_volume("TRIANGLE_TRI7_1", :Tri7), 1/2)
+    @test isapprox(calculate_volume("SQUARE_QUAD4_1", :Quad4), 2^2)
+    @test isapprox(calculate_volume("SQUARE_QUAD8_1", :Quad8), 2^2)
+    @test isapprox(calculate_volume("SQUARE_QUAD9_1", :Quad9), 2^2)
+    @test isapprox(calculate_volume("TETRA_TET4_1", :Tet4), 1/6)
+    @test isapprox(calculate_volume("TETRA_TET10_1", :Tet10), 1/6)
+#   @test isapprox(calculate_volume("TETRA_TET14_1", :Tet14), 1/6)
+    @test isapprox(calculate_volume("CUBE_HEX8_1", :Hex8), 2^3)
+    @test isapprox(calculate_volume("CUBE_HEX20_1", :Hex20), 2^3)
+    @test isapprox(calculate_volume("CUBE_HEX27_1", :Hex27), 2^3)
+    @test isapprox(calculate_volume("WEDGE_WEDGE6_1", :Wedge6), 1)
+#   @test isapprox(calculate_volume("WEDGE_WEDGE15_1", :Wedge15, 1/2))
+#   @test isapprox(calculate_volume("PYRAMID_PYRAMID5_1", :Pyramid5, ?))
+#   @test isapprox(calculate_volume("PYRAMID_PYRAMID13_1", :Pyramid13, ?))
+end
+
+@testset "read nodal field from code aster result file" begin
+    rmedfile = joinpath(datadir, "rings.rmed")
+    rmed = JuliaFEM.Preprocess.RMEDFile(rmedfile)
+    temp = JuliaFEM.Preprocess.aster_read_data(rmed, "TEMP")
+    @test isapprox(temp[15], 1.0)
+    @test isapprox(temp[95], 2.0)
+end
diff --git a/test/test_preprocess/block_2d_1elem_quad4.med b/test/test_preprocess/block_2d_1elem_quad4.med
new file mode 100644
index 0000000..7948982
Binary files /dev/null and b/test/test_preprocess/block_2d_1elem_quad4.med differ
diff --git a/test/test_preprocess/primitives.med b/test/test_preprocess/primitives.med
new file mode 100644
index 0000000..8393351
Binary files /dev/null and b/test/test_preprocess/primitives.med differ
diff --git a/test/test_preprocess/rings.rmed b/test/test_preprocess/rings.rmed
new file mode 100644
index 0000000..e0f1a0d
Binary files /dev/null and b/test/test_preprocess/rings.rmed differ
diff --git a/test/test_preprocess_aster_reader.jl b/test/test_preprocess_aster_reader.jl
deleted file mode 100644
index ba01619..0000000
--- a/test/test_preprocess_aster_reader.jl
+++ /dev/null
@@ -1,158 +0,0 @@
-# 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.Testing
-
-@testset "read ascii mesh" begin
-    mesh = """
-    COOR_2D
-    N1          0.0 0.0
-    N2          1.0 0.0
-    N3          1.0 1.0
-    N4          0.0 1.0
-    FINSF
-
-    QUAD4
-    E1         N1 N2 N3 N4
-    FINSF
-
-    SEG2
-    E2         N3 N4
-    FINSF
-
-    GROUP_NO NOM=NALL
-    N1 N2
-    FINSF
-
-    GROUP_MA NOM=BODY1
-    E1 E2
-    FINSF
-
-    FIN
-    """
-
-#   m = parse(mesh, Val{:CODE_ASTER_MAIL})
-#   @test m["nodes"]["N1"] == [0.0, 0.0]
-#   @test m["elements"]["E1"] == ["QUAD4", ["N1", "N2", "N3", "N4"]]
-#   @test m["elements"]["E2"] == ["SEG2", ["N3", "N4"]]
-#   @test m["elsets"]["BODY1"] == ["E1", "E2"]
-#   @test m["nsets"]["NALL"] == ["N1", "N2"]
-end
-
-@testset "parse nodes" begin
-    section = """
-    N9  2.0 3.0 4.0
-    COOR_3D
-    N1          0.0 0.0 0.0
-    N2          1.0 0.0 0.0
-    N3          1.0 1.0 0.0
-    N4          0.0 1.0 0.0
-    N5          0.0 0.0 1.0
-    N6          1.0 0.0 1.0
-    N7          1.0 1.0 1.0
-    N8          0.0 1.0 1.0
-    FINSF
-    absdflasdf
-    N12 3.0 4.0 5.0 6.0
-    N13 3.0 4.0 5.0
-    """
-    nodes = aster_parse_nodes(section)
-    @test nodes[1] == Float64[0.0, 0.0, 0.0]
-    @test nodes[8] == Float64[0.0, 1.0, 1.0]
-    @test length(nodes) == 8
-end
-
-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)
-    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
-end
-
-function calculate_volume(mesh_name, eltype)
-    fn = Pkg.dir("JuliaFEM") * "/test/testdata/primitives.med"
-    mesh = aster_read_mesh(fn, mesh_name)
-    elements = create_elements(mesh; element_type=eltype)
-    V = 0.0
-    time = 0.0
-    for element in elements
-        for ip in get_integration_points(element)
-            detJ = element(ip, time, Val{:detJ})
-            detJ > 0 || warn("negative determinant for element $eltype !")
-            V += ip.weight*detJ
-        end
-    end
-    info("volume of $eltype is $V")
-    return V
-end
-
-@testset "calculate volume for 1 element models" begin
-    @test isapprox(calculate_volume("TRIANGLE_TRI3_1", :Tri3), 1/2)
-    @test isapprox(calculate_volume("TRIANGLE_TRI6_1", :Tri6), 1/2)
-    @test isapprox(calculate_volume("TRIANGLE_TRI7_1", :Tri7), 1/2)
-    @test isapprox(calculate_volume("SQUARE_QUAD4_1", :Quad4), 2^2)
-    @test isapprox(calculate_volume("SQUARE_QUAD8_1", :Quad8), 2^2)
-    @test isapprox(calculate_volume("SQUARE_QUAD9_1", :Quad9), 2^2)
-    @test isapprox(calculate_volume("TETRA_TET4_1", :Tet4), 1/6)
-    @test isapprox(calculate_volume("TETRA_TET10_1", :Tet10), 1/6)
-#   @test isapprox(calculate_volume("TETRA_TET14_1", :Tet14), 1/6)
-    @test isapprox(calculate_volume("CUBE_HEX8_1", :Hex8), 2^3)
-    @test isapprox(calculate_volume("CUBE_HEX20_1", :Hex20), 2^3)
-    @test isapprox(calculate_volume("CUBE_HEX27_1", :Hex27), 2^3)
-    @test isapprox(calculate_volume("WEDGE_WEDGE6_1", :Wedge6), 1)
-#   @test isapprox(calculate_volume("WEDGE_WEDGE15_1", :Wedge15, 1/2))
-#   @test isapprox(calculate_volume("PYRAMID_PYRAMID5_1", :Pyramid5, ?))
-#   @test isapprox(calculate_volume("PYRAMID_PYRAMID13_1", :Pyramid13, ?))
-end
-
-@testset "get nodal field from aster file" begin
-    fn = Pkg.dir("JuliaFEM") * "/test/testdata/rings.rmed"
-    medfile = JuliaFEM.Preprocess.RMEDFile(fn)
-    temp = JuliaFEM.Preprocess.aster_read_data(medfile, "TEMP")
-    info("temp = $temp")
-    # more like functional testing, results are what they are,
-    # we're happy to just have some results
-    @test true
-end
diff --git a/test/test_problems_contact_2d.jl b/test/test_problems_contact_2d.jl
new file mode 100644
index 0000000..ccd973d
--- /dev/null
+++ b/test/test_problems_contact_2d.jl
@@ -0,0 +1,257 @@
+# 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.Testing
+
+datadir = first(splitext(basename(@__FILE__)))
+
+# from fenet d3613 advanced finite element contact benchmarks
+# a = 6.21 mm, pmax = 3585 MPa
+# this is a very sparse mesh and for that reason pmax is not very accurate
+# (only 6 elements in -20 .. 20 mm contact zone, 3 elements in contact
+@testset "hertz contact, full 2d model, linear elements, curved slave surface" begin
+    meshfile = joinpath(datadir, "hertz_2d_full.med")
+    mesh = aster_read_mesh(meshfile)
+
+    upper = Problem(Elasticity, "CYLINDER", 2)
+    upper.properties.formulation = :plane_strain
+    upper.elements = create_elements(mesh, "CYLINDER")
+    update!(upper, "youngs modulus", 70.0e3)
+    update!(upper, "poissons ratio", 0.3)
+
+    lower = Problem(Elasticity, "BLOCK", 2)
+    lower.properties.formulation = :plane_strain
+    lower.elements = create_elements(mesh, "BLOCK")
+    update!(lower, "youngs modulus", 210.0e3)
+    update!(lower, "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, "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, "displacement 1", 0.0)
+
+    nid = find_nearest_node(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, "displacement traction force 2", -35.0e3)
+
+    contact = Problem(Contact, "contact between block and cylinder", 2, "displacement")
+    contact.properties.rotate_normals = true
+    contact.properties.finite_sliding = false
+    contact.properties.friction = false
+    contact.properties.use_forwarddiff = false
+    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)
+    solver()
+    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)
+    # 12 % error in maximum pressure
+    # 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
+    info("2d hertz: Rn = $Rn, Rt = $Rt")
+    info("2d hertz: maximum pressure pmax = ", maximum(pres))
+    @test isapprox(maximum(pres), 3585.0; rtol = 0.13)
+    # under 0.15 % error in reaction force
+    @test isapprox(Rn, 35.0e3; rtol=0.020)
+    @test isapprox(Rt, 0.0; atol=200.0)
+end
+
+@testset "hertz contact, full 2d model, linear elements, flat slave surface" begin
+    meshfile = joinpath(datadir, "hertz_2d_full.med")
+    mesh = aster_read_mesh(meshfile)
+
+    upper = Problem(Elasticity, "CYLINDER", 2)
+    upper.properties.formulation = :plane_strain
+    upper.elements = create_elements(mesh, "CYLINDER")
+    update!(upper, "youngs modulus", 70.0e3)
+    update!(upper, "poissons ratio", 0.3)
+
+    lower = Problem(Elasticity, "BLOCK", 2)
+    lower.properties.formulation = :plane_strain
+    lower.elements = create_elements(mesh, "BLOCK")
+    update!(lower, "youngs modulus", 210.0e3)
+    update!(lower, "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, "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, "displacement 1", 0.0)
+
+    nid = find_nearest_node(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, "displacement traction force 2", -35.0e3)
+
+    contact = Problem(Contact, "contact between block and cylinder", 2, "displacement")
+    contact.properties.rotate_normals = true
+    contact.properties.finite_sliding = false
+    contact.properties.friction = false
+    contact.properties.use_forwarddiff = false
+    contact_slave_elements = create_elements(mesh, "BLOCK_TO_CYLINDER")
+    contact_master_elements = create_elements(mesh, "CYLINDER_TO_BLOCK")
+    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)
+    solver()
+    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)
+    # 12 % error in maximum pressure
+    # 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
+    info("2d hertz: Rn = $Rn, Rt = $Rt")
+    info("2d hertz: maximum pressure pmax = ", maximum(pres))
+    @test isapprox(maximum(pres), 3585.0; rtol = 0.13)
+    # under 0.15 % error in reaction force
+    @test isapprox(Rn, 35.0e3; rtol=0.020)
+    @test isapprox(Rt, 0.0; atol=200.0)
+end
+
+function get_model()
+    meshfile = joinpath(datadir, "block_2d.med")
+    mesh = aster_read_mesh(meshfile)
+    println(mesh.nodes[1])
+
+    upper = Problem(mesh, Elasticity, "UPPER", 2)
+    lower = Problem(mesh, Elasticity, "LOWER", 2)
+
+    for body in [upper, lower]
+        body.properties.formulation = :plane_stress
+        update!(body, "youngs modulus", 288.0)
+        update!(body, "poissons ratio", 1/3)
+    end
+
+    load = Problem(mesh, Elasticity, "UPPER_TOP", 2)
+    load.properties.formulation = :plane_stress
+    update!(load, "displacement traction force 2", -28.8)
+    bc1 = Problem(mesh, Dirichlet, "LOWER_BOTTOM", 2, "displacement")
+    update!(bc1, "displacement 2", 0.0)
+    bc2 = Problem(mesh, Dirichlet, "LOWER_LEFT", 2, "displacement")
+    update!(bc2, "displacement 1", 0.0)
+    bc3 = Problem(mesh, Dirichlet, "UPPER_LEFT", 2, "displacement")
+    update!(bc3, "displacement 1", 0.0)
+
+    interface = Problem(Contact, "interface", 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]
+    interface.properties.rotate_normals = true
+    
+    # in LOWER_LEFT we have node belonging also to contact interface
+    # let's remove it from dirichlet bc
+    create_node_set_from_element_set!(mesh, "LOWER_LEFT")
+    nid = find_nearest_node(mesh, [0.0, 0.5]; node_set="LOWER_LEFT")
+    coords = mesh.nodes[nid]
+    info("nearest node to (0.0, 0.5) = $nid, coordinates = $coords")
+    dofs = [2*(nid-1)+1, 2*(nid-1)+2]
+    info("removing nid $nid, dofs $dofs from LOWER_LEFT")
+    push!(bc2.assembly.removed_dofs, dofs...)
+
+    solver = Solver(Nonlinear)
+    #push!(solver, upper, lower, load, bc1, interface)
+    push!(solver, upper, lower, load, bc1, bc2, bc3, interface)
+    return solver
+end
+
+@testset "small sliding 2d patch test, linear Seg2 elements, standard basis" begin
+
+    solver = get_model()
+    interface = solver["interface"]
+    solver()
+
+    node_ids, displacement = get_nodal_vector(interface.elements, "displacement", 0.0)
+    node_ids, geometry = get_nodal_vector(interface.elements, "geometry", 0.0)
+    node_ids, reaction_force = get_nodal_vector(get_slave_elements(interface), "reaction force", 0.0)
+    u2 = [u[2] for u in displacement]
+    f2 = [f[2] for f in reaction_force]
+    maxabsu2 = maximum(abs(u2))
+    stdabsu2 = std(abs(u2))
+    info("max(abs(u2)) = $maxabsu2, std(abs(u2)) = $stdabsu2")
+    @test isapprox(stdabsu2, 0.0; atol=1.0e-12)
+    maxabsf2 = maximum(abs(f2))
+    stdabsf2 = std(abs(f2))
+    info("max(abs(f2)) = $maxabsf2, std(abs(f2)) = $stdabsf2")
+    @test isapprox(stdabsf2, 0.0; atol=1.0e-12)
+    @test isapprox(mean(abs(f2)), 28.8; atol=1.0e-12)
+end
+
+@testset "small sliding 2d patch test, linear Seg2 elements, dual basis" begin
+
+    solver = get_model()
+    interface = solver["interface"]
+    interface.properties.dual_basis = true
+    solver()
+
+    node_ids, displacement = get_nodal_vector(interface.elements, "displacement", 0.0)
+    node_ids, geometry = get_nodal_vector(interface.elements, "geometry", 0.0)
+    node_ids, reaction_force = get_nodal_vector(get_slave_elements(interface), "reaction force", 0.0)
+    u2 = [u[2] for u in displacement]
+    f2 = [f[2] for f in reaction_force]
+    maxabsu2 = maximum(abs(u2))
+    stdabsu2 = std(abs(u2))
+    info("max(abs(u2)) = $maxabsu2, std(abs(u2)) = $stdabsu2")
+    @test isapprox(stdabsu2, 0.0; atol=1.0e-12)
+    maxabsf2 = maximum(abs(f2))
+    stdabsf2 = std(abs(f2))
+    info("max(abs(f2)) = $maxabsf2, std(abs(f2)) = $stdabsf2")
+    @test isapprox(stdabsf2, 0.0; atol=1.0e-12)
+    @test isapprox(mean(abs(f2)), 28.8; atol=1.0e-12)
+end
diff --git a/test/test_problems_contact_2d/block_2d.med b/test/test_problems_contact_2d/block_2d.med
new file mode 100644
index 0000000..55e8b7a
Binary files /dev/null and b/test/test_problems_contact_2d/block_2d.med differ
diff --git a/test/testdata/hertz_2d_full.med b/test/test_problems_contact_2d/hertz_2d_full.med
similarity index 100%
rename from test/testdata/hertz_2d_full.med
rename to test/test_problems_contact_2d/hertz_2d_full.med
diff --git a/test/test_problems_contact_2d_autodiff.jl b/test/test_problems_contact_2d_autodiff.jl
new file mode 100644
index 0000000..6f26896
--- /dev/null
+++ b/test/test_problems_contact_2d_autodiff.jl
@@ -0,0 +1,92 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
+
+using JuliaFEM
+using JuliaFEM.Preprocess
+using JuliaFEM.Postprocess
+using JuliaFEM.Testing
+
+datadir = first(splitext(basename(@__FILE__)))
+
+function get_model()
+    meshfile = joinpath(datadir, "block_2d.med")
+    mesh = aster_read_mesh(meshfile)
+    println(mesh.nodes[1])
+
+    upper = Problem(mesh, Elasticity, "UPPER", 2)
+    lower = Problem(mesh, Elasticity, "LOWER", 2)
+
+    for body in [upper, lower]
+        body.properties.formulation = :plane_stress
+        update!(body, "youngs modulus", 288.0)
+        update!(body, "poissons ratio", 1/3)
+    end
+
+    load = Problem(mesh, Elasticity, "UPPER_TOP", 2)
+    load.properties.formulation = :plane_stress
+    update!(load, "displacement traction force 2", 0.0 => 0.0)
+    update!(load, "displacement traction force 2", 1.0 => -28.8)
+    bc1 = Problem(mesh, Dirichlet, "LOWER_BOTTOM", 2, "displacement")
+    update!(bc1, "displacement 2", 0.0)
+    bc2 = Problem(mesh, Dirichlet, "LOWER_LEFT", 2, "displacement")
+    update!(bc2, "displacement 1", 0.0)
+    bc3 = Problem(mesh, Dirichlet, "UPPER_LEFT", 2, "displacement")
+    update!(bc3, "displacement 1", 0.0)
+
+    interface = Problem(Contact, "interface", 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]
+    interface.properties.rotate_normals = true
+    
+    # in LOWER_LEFT we have node belonging also to contact interface
+    # let's remove it from dirichlet bc
+    create_node_set_from_element_set!(mesh, "LOWER_LEFT")
+    nid = find_nearest_node(mesh, [0.0, 0.5]; node_set="LOWER_LEFT")
+    coords = mesh.nodes[nid]
+    info("nearest node to (0.0, 0.5) = $nid, coordinates = $coords")
+    dofs = [2*(nid-1)+1, 2*(nid-1)+2]
+    info("removing nid $nid, dofs $dofs from LOWER_LEFT")
+    push!(bc2.assembly.removed_dofs, dofs...)
+
+    solver = Solver(Nonlinear)
+    push!(solver, upper, lower, load, bc1, bc2, bc3, interface)
+    return solver
+end
+
+@testset "finite sliding 2d patch test, linear Seg2 elements, standard basis" begin
+
+    solver = get_model()
+    interface = solver["interface"]
+    upper = solver["UPPER"]
+    lower = solver["LOWER"]
+    for body in [upper, lower]
+        body.properties.geometric_stiffness = true
+        body.properties.finite_strain = true
+    end
+    interface.properties.finite_sliding = true
+    interface.properties.use_forwarddiff = true
+
+    for time in [0.0, 1/3, 2/3, 1.0]
+        interface.properties.iteration = 1
+        solver.time = time
+        solver()
+    end
+
+    node_ids, displacement = get_nodal_vector(interface.elements, "displacement", 1.0)
+    node_ids, geometry = get_nodal_vector(interface.elements, "geometry", 1.0)
+    node_ids, reaction_force = get_nodal_vector(get_slave_elements(interface), "reaction force", 1.0)
+    u2 = [u[2] for u in displacement]
+    f2 = [f[2] for f in reaction_force]
+    maxabsu2 = maximum(abs(u2))
+    stdabsu2 = std(abs(u2))
+    info("max(abs(u2)) = $maxabsu2, std(abs(u2)) = $stdabsu2")
+    @test isapprox(stdabsu2, 0.0; atol=1.0e-12)
+    maxabsf2 = maximum(abs(f2))
+    stdabsf2 = std(abs(f2))
+    info("max(abs(f2)) = $maxabsf2, std(abs(f2)) = $stdabsf2")
+    @test isapprox(stdabsf2, 0.0; atol=1.0e-12)
+    # for linear case pressure 28.8
+    @test isapprox(mean(abs(f2)), 27.76616800689944; rtol=1.0e-3)
+end
diff --git a/test/test_problems_contact_2d_autodiff/block_2d.med b/test/test_problems_contact_2d_autodiff/block_2d.med
new file mode 100644
index 0000000..55e8b7a
Binary files /dev/null and b/test/test_problems_contact_2d_autodiff/block_2d.med differ
diff --git a/test/test_problems_contact_3d.jl b/test/test_problems_contact_3d.jl
new file mode 100644
index 0000000..2c27f9d
--- /dev/null
+++ b/test/test_problems_contact_3d.jl
@@ -0,0 +1,131 @@
+# 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.Testing
+using JuliaFEM.Abaqus: create_surface_elements
+
+tet4_meshfile = "test_problems_contact_3d/tet4.inp"
+tet10_meshfile = "test_problems_contact_3d/tet10.inp"
+
+function get_model(meshfile)
+    mesh = abaqus_read_mesh(meshfile)
+    
+    upper = Problem(Elasticity, "UPPER", 3)
+    upper.elements = create_elements(mesh, "UPPER")
+    update!(upper, "youngs modulus", 3*288.0)
+    update!(upper, "poissons ratio", 1/3)
+
+    lower = Problem(Elasticity, "LOWER", 3)
+    lower.elements = create_elements(mesh, "LOWER")
+    update!(lower, "youngs modulus", 288.0)
+    update!(lower, "poissons ratio", 1/3)
+
+    bc_upper = Problem(Dirichlet, "UPPER_TOP", 3, "displacement")
+    bc_upper.elements = create_surface_elements(mesh, "UPPER_TOP")
+    update!(bc_upper, "displacement 3", -0.4)
+
+    bc_lower = Problem(Dirichlet, "LOWER_BOTTOM", 3, "displacement")
+    bc_lower.elements = create_surface_elements(mesh, "LOWER_BOTTOM")
+    update!(bc_lower, "displacement 3", 0.0)
+    
+    # point-wise boundary conditions to prevent free body move
+    nid1 = find_nearest_nodes(mesh, [0.0, 0.0, 0.0])[1]
+    nid2 = find_nearest_nodes(mesh, [1.0, 0.0, 0.0])[1]
+    nid3 = find_nearest_nodes(mesh, [0.0, 1.0, 0.0])[1]
+    nid4 = find_nearest_nodes(mesh, [0.0, 0.0, 1.0])[1]
+    nid5 = find_nearest_nodes(mesh, [1.0, 0.0, 1.0])[1]
+    nid6 = find_nearest_nodes(mesh, [0.0, 1.0, 1.0])[1]
+
+    bc_sym13 = Problem(Dirichlet, "SYM13", 3, "displacement")
+    # nodes in X2=0 plane
+    bc_sym13.elements = [Element(Poi1, [j]) for j in [nid1, nid2, nid4, nid5]]
+    update!(bc_sym13, "geometry", mesh.nodes)
+    update!(bc_sym13, "displacement 2", 0.0)
+
+    bc_sym23 = Problem(Dirichlet, "SYM23", 3, "displacement")
+    # nodes in X1=0 plane
+    bc_sym23.elements = [Element(Poi1, [j]) for j in [nid1, nid3, nid4, nid6]]
+    update!(bc_sym23, "geometry", mesh.nodes)
+    update!(bc_sym23, "displacement 1", 0.0)
+
+    interface = Problem(Contact, "LOWER_TO_UPPER", 3, "displacement")
+    interface_slave_elements = create_surface_elements(mesh, "LOWER_TO_UPPER")
+    interface_master_elements = create_surface_elements(mesh, "UPPER_TO_LOWER")
+    update!(interface_slave_elements, "master elements", interface_master_elements)
+    interface.elements = [interface_slave_elements; interface_master_elements]
+    interface.properties.contact_state_in_first_iteration = :AUTO
+    
+    #append!(bc_sym13.assembly.removed_dofs, [1316, 1319, 1388, 1358])
+    #append!(bc_sym23.assembly.removed_dofs, [1492, 1627, 1387, 1597])
+    #append!(interface.assembly.removed_dofs, [1316, 1319, 1358, 1387, 1388, 1492, 1597, 1627])
+
+    solver = NonlinearSolver(upper, lower, bc_upper, bc_lower, bc_sym13, bc_sym23, interface)
+#    solver.properties.max_iterations = 5
+
+    return solver
+end
+
+@testset "small sliding contact patch test, tet4 + standard basis" begin
+    solver = get_model(tet4_meshfile)
+    solver.xdmf = Xdmf("contact_sl_lin_disp_results"; overwrite=true)
+    interface = solver["LOWER_TO_UPPER"]
+    interface.properties.dual_basis = false
+    solver()
+    node_ids, displacement = get_nodal_vector(interface.elements, "displacement", 0.0)
+    node_ids, geometry = get_nodal_vector(interface.elements, "geometry", 0.0)
+    # node_ids, pressure = get_nodal_vector(interface.elements, "contact pressure", 0.0)
+    u3 = [u[3] for u in displacement]
+    maxabsu3 = maximum(abs(u3))
+    stdabsu3 = std(abs(u3))
+    info("max(abs(u3)) = $maxabsu3, std(abs(u3)) = $stdabsu3")
+    @test isapprox(stdabsu3, 0.0; atol=1.0e-12)
+end
+
+@testset "small sliding contact patch test, tet4 + dual basis" begin
+    solver = get_model(tet4_meshfile)
+    solver.xdmf = Xdmf("contact_dl_lin_disp_results"; overwrite=true)
+    interface = solver["LOWER_TO_UPPER"]
+    interface.properties.dual_basis = true
+    solver()
+    node_ids, displacement = get_nodal_vector(interface.elements, "displacement", 0.0)
+    node_ids, geometry = get_nodal_vector(interface.elements, "geometry", 0.0)
+    u3 = [u[3] for u in displacement]
+    maxabsu3 = maximum(abs(u3))
+    stdabsu3 = std(abs(u3))
+    info("max(abs(u3)) = $maxabsu3, std(abs(u3)) = $stdabsu3")
+    @test isapprox(stdabsu3, 0.0; atol=1.0e-12)
+end
+
+@testset "small sliding contact patch test, tet10 + standard basis" begin
+    solver = get_model(tet10_meshfile)
+    solver.xdmf = Xdmf("contact_sl_quad_disp_results"; overwrite=true)
+    interface = solver["LOWER_TO_UPPER"]
+    interface.properties.dual_basis = false
+    solver()
+    node_ids, displacement = get_nodal_vector(interface.elements, "displacement", 0.0)
+    node_ids, geometry = get_nodal_vector(interface.elements, "geometry", 0.0)
+    u3 = [u[3] for u in displacement]
+    maxabsu3 = maximum(abs(u3))
+    stdabsu3 = std(abs(u3))
+    info("max(abs(u3)) = $maxabsu3, std(abs(u3)) = $stdabsu3")
+    @test isapprox(stdabsu3, 0.0; atol=1.0e-10)
+end
+
+@testset "small sliding contact patch test, tet10 + dual basis, alpha=0.2" begin
+    solver = get_model(tet10_meshfile)
+    solver.xdmf = Xdmf("contact_sl_quad_disp_results"; overwrite=true)
+    interface = solver["LOWER_TO_UPPER"]
+    interface.properties.dual_basis = true
+    interface.properties.alpha = 0.2
+    solver()
+    node_ids, displacement = get_nodal_vector(interface.elements, "displacement", 0.0)
+    node_ids, geometry = get_nodal_vector(interface.elements, "geometry", 0.0)
+    u3 = [u[3] for u in displacement]
+    maxabsu3 = maximum(abs(u3))
+    stdabsu3 = std(abs(u3))
+    info("max(abs(u3)) = $maxabsu3, std(abs(u3)) = $stdabsu3")
+    @test isapprox(stdabsu3, 0.0; atol=1.0e-10)
+end
diff --git a/test/test_problems_contact_3d/tet10.inp b/test/test_problems_contact_3d/tet10.inp
new file mode 100644
index 0000000..6579ac9
--- /dev/null
+++ b/test/test_problems_contact_3d/tet10.inp
@@ -0,0 +1,1354 @@
+**NSET COUNT = 607
+*NODE
+1,                0.74534,        0.75046,        0.82163
+2,                0.62267,        0.75023,        0.91082
+3,                0.62267,        0.87523,        0.91082
+4,                0.74767,        0.75023,        0.91082
+5,                0.87500,        0.50000,        0.87500
+6,                0.74767,        0.62523,        0.91082
+7,                0.87267,        0.62523,        0.91082
+8,                0.87267,        0.62523,        0.78582
+9,                0.87267,        0.75023,        0.91082
+10,               0.87267,        0.75023,        0.78582
+11,               0.87500,        0.62500,        0.62500
+12,               0.87500,        0.75000,        0.62500
+13,               0.74767,        0.75023,        0.66082
+14,               0.66838,        0.54400,        0.72860
+15,               0.87500,        0.50000,        0.62500
+16,               0.70919,        0.52200,        0.61430
+17,               0.83419,        0.52200,        0.73930
+18,               0.70919,        0.64700,        0.61430
+19,               0.62500,        0.87500,        0.62500
+20,               0.50000,        0.87500,        0.62500
+21,               0.62267,        0.75023,        0.66082
+22,               0.62267,        0.87523,        0.78582
+23,               0.75000,        0.87500,        0.62500
+24,               0.74767,        0.87523,        0.78582
+25,               0.87500,        0.87500,        0.75000
+26,               0.87500,        0.87500,        0.87500
+27,               0.74767,        0.87523,        0.91082
+28,               0.58419,        0.52200,        0.86430
+29,               0.58419,        0.64700,        0.86430
+30,               0.70919,        0.52200,        0.86430
+31,               0.58419,        0.39700,        0.86430
+32,               0.70919,        0.39700,        0.86430
+33,               0.87500,        0.37500,        0.87500
+34,               0.87500,        0.25000,        0.87500
+35,               0.83419,        0.39700,        0.73930
+36,               0.58419,        0.39700,        0.61430
+37,               0.70919,        0.39700,        0.61430
+38,               0.58419,        0.52200,        0.61430
+39,               0.87500,        0.25000,        0.62500
+40,               0.87500,        0.37500,        0.62500
+41,               0.70686,        0.64723,        0.77512
+42,               0.87500,        0.87500,        0.62500
+43,               0.58419,        0.64700,        0.61430
+44,               0.75000,        0.12500,        0.62500
+45,               0.87500,        0.12500,        0.62500
+46,               0.87500,        0.12500,        0.87500
+47,               0.87500,        0.12500,        0.75000
+48,               0.37836,        0.85595,        0.74983
+49,               0.43918,        0.80298,        0.87491
+50,               0.56185,        0.80321,        0.78573
+51,               0.52337,        0.69998,        0.73922
+52,               0.43918,        0.67798,        0.87491
+53,               0.31418,        0.80298,        0.87491
+54,               0.31418,        0.67798,        0.87491
+55,               0.12500,        0.62500,        0.87500
+56,               0.12500,        0.75000,        0.87500
+57,               0.18918,        0.80298,        0.74991
+58,               0.12500,        0.50000,        0.87500
+59,               0.18918,        0.67798,        0.74991
+60,               0.43918,        0.67798,        0.62491
+61,               0.43918,        0.80298,        0.62491
+62,               0.43918,        0.92798,        0.74991
+63,               0.43918,        0.92798,        0.62491
+64,               0.31418,        0.92798,        0.62491
+65,               0.31418,        0.80298,        0.62491
+66,               0.31418,        0.67798,        0.62491
+67,               0.12500,        0.62500,        0.62500
+68,               0.12500,        0.50000,        0.62500
+69,               0.12500,        0.75000,        0.62500
+70,               0.12500,        0.87500,        0.62500
+71,               0.31418,        0.92798,        0.74991
+72,               0.31418,        0.92798,        0.87491
+73,               0.25000,        0.87500,        0.87500
+74,               0.43918,        0.92798,        0.87491
+75,               0.12500,        0.87500,        0.75000
+76,               0.62345,        0.16661,        0.75334
+77,               0.50000,        0.12500,        0.87500
+78,               0.56172,        0.08330,        0.87667
+79,               0.56172,        0.08330,        0.75167
+80,               0.56172,        0.20830,        0.87667
+81,               0.68672,        0.08330,        0.75167
+82,               0.68672,        0.08330,        0.87667
+83,               0.68672,        0.20830,        0.87667
+84,               0.64591,        0.35531,        0.74097
+85,               0.81172,        0.20830,        0.75167
+86,               0.68672,        0.20830,        0.62667
+87,               0.56172,        0.20830,        0.62667
+88,               0.68672,        0.08330,        0.62667
+89,               0.56172,        0.08330,        0.62667
+90,               0.22841,        0.22025,        0.65164
+91,               0.23920,        0.23513,        0.57582
+92,               0.36420,        0.11013,        0.57582
+93,               0.36420,        0.23513,        0.57582
+94,               0.41565,        0.41697,        0.69068
+95,               0.33283,        0.45849,        0.84534
+96,               0.39700,        0.63646,        0.72025
+97,               0.45783,        0.45849,        0.84534
+98,               0.20783,        0.45849,        0.72034
+99,               0.33283,        0.33349,        0.84534
+100,              0.45783,        0.33349,        0.84534
+101,              0.33283,        0.33349,        0.59534
+102,              0.45783,        0.33349,        0.59534
+103,              0.33283,        0.45849,        0.59534
+104,              0.51955,        0.29179,        0.72201
+105,              0.54202,        0.48049,        0.70964
+106,              0.45783,        0.45849,        0.59534
+107,              0.37500,        0.12500,        0.87500
+108,              0.25000,        0.12500,        0.87500
+109,              0.23920,        0.23513,        0.82582
+110,              0.36420,        0.11013,        0.70082
+111,              0.23920,        0.11013,        0.70082
+112,              0.12500,        0.12500,        0.75000
+113,              0.12500,        0.25000,        0.87500
+114,              0.11420,        0.23513,        0.70082
+115,              0.12500,        0.37500,        0.87500
+116,              0.11420,        0.36013,        0.70082
+117,              0.32203,        0.31861,        0.67116
+118,              0.36420,        0.23513,        0.82582
+119,              0.42593,        0.19343,        0.70249
+120,              0.12500,        0.12500,        0.87500
+121,              0.23920,        0.11013,        0.57582
+122,              0.12500,        0.12500,        0.62500
+123,              0.11420,        0.23513,        0.57582
+124,              0.11420,        0.36013,        0.57582
+125,              0.23920,        0.36013,        0.57582
+126,              0.12500,        0.87500,        0.87500
+127,              1.00000,        0.00000,        0.75000
+128,              1.00000,        0.00000,        1.00000
+129,              0.75000,        0.00000,        1.00000
+130,              1.00000,        0.00000,        0.87500
+131,              0.87500,        0.00000,        1.00000
+132,              0.87500,        0.00000,        0.87500
+133,              0.75000,        0.00000,        0.75000
+134,              0.75000,        0.00000,        0.87500
+135,              0.87500,        0.00000,        0.75000
+136,              1.00000,        0.00000,        0.50000
+137,              1.00000,        0.00000,        0.62500
+138,              0.87500,        0.00000,        0.62500
+139,              0.75000,        0.00000,        0.50000
+140,              0.75000,        0.00000,        0.62500
+141,              0.87500,        0.00000,        0.50000
+142,              0.50000,        0.00000,        1.00000
+143,              0.62500,        0.00000,        1.00000
+144,              0.62500,        0.00000,        0.87500
+145,              0.50000,        0.00000,        0.75000
+146,              0.50000,        0.00000,        0.87500
+147,              0.62500,        0.00000,        0.75000
+148,              0.62500,        0.00000,        0.62500
+149,              0.50000,        0.00000,        0.50000
+150,              0.50000,        0.00000,        0.62500
+151,              0.62500,        0.00000,        0.50000
+152,              0.25000,        0.00000,        1.00000
+153,              0.37500,        0.00000,        1.00000
+154,              0.37500,        0.00000,        0.87500
+155,              0.25000,        0.00000,        0.75000
+156,              0.25000,        0.00000,        0.87500
+157,              0.37500,        0.00000,        0.75000
+158,              0.37500,        0.00000,        0.62500
+159,              0.25000,        0.00000,        0.50000
+160,              0.25000,        0.00000,        0.62500
+161,              0.37500,        0.00000,        0.50000
+162,              0.00000,        0.00000,        1.00000
+163,              0.12500,        0.00000,        1.00000
+164,              0.12500,        0.00000,        0.87500
+165,              0.00000,        0.00000,        0.75000
+166,              0.00000,        0.00000,        0.87500
+167,              0.12500,        0.00000,        0.75000
+168,              0.12500,        0.00000,        0.62500
+169,              0.00000,        0.00000,        0.50000
+170,              0.00000,        0.00000,        0.62500
+171,              0.12500,        0.00000,        0.50000
+172,              0.75000,        1.00000,        1.00000
+173,              1.00000,        1.00000,        1.00000
+174,              1.00000,        1.00000,        0.75000
+175,              0.87500,        1.00000,        1.00000
+176,              1.00000,        1.00000,        0.87500
+177,              0.87500,        1.00000,        0.87500
+178,              0.75000,        1.00000,        0.75000
+179,              0.75000,        1.00000,        0.87500
+180,              0.87500,        1.00000,        0.75000
+181,              1.00000,        1.00000,        0.50000
+182,              1.00000,        1.00000,        0.62500
+183,              0.87500,        1.00000,        0.62500
+184,              0.75000,        1.00000,        0.50000
+185,              0.75000,        1.00000,        0.62500
+186,              0.87500,        1.00000,        0.50000
+187,              0.50000,        1.00000,        1.00000
+188,              0.62500,        1.00000,        1.00000
+189,              0.62500,        1.00000,        0.87500
+190,              0.50000,        1.00000,        0.75000
+191,              0.50000,        1.00000,        0.87500
+192,              0.62500,        1.00000,        0.75000
+193,              0.62500,        1.00000,        0.62500
+194,              0.50000,        1.00000,        0.50000
+195,              0.50000,        1.00000,        0.62500
+196,              0.62500,        1.00000,        0.50000
+197,              0.25000,        1.00000,        1.00000
+198,              0.37500,        1.00000,        1.00000
+199,              0.37500,        1.00000,        0.87500
+200,              0.25000,        1.00000,        0.75000
+201,              0.25000,        1.00000,        0.87500
+202,              0.37500,        1.00000,        0.75000
+203,              0.37500,        1.00000,        0.62500
+204,              0.25000,        1.00000,        0.50000
+205,              0.25000,        1.00000,        0.62500
+206,              0.37500,        1.00000,        0.50000
+207,              0.00000,        1.00000,        1.00000
+208,              0.12500,        1.00000,        1.00000
+209,              0.12500,        1.00000,        0.87500
+210,              0.00000,        1.00000,        0.75000
+211,              0.00000,        1.00000,        0.87500
+212,              0.12500,        1.00000,        0.75000
+213,              0.12500,        1.00000,        0.62500
+214,              0.00000,        1.00000,        0.50000
+215,              0.00000,        1.00000,        0.62500
+216,              0.12500,        1.00000,        0.50000
+217,              1.00000,        0.75000,        1.00000
+218,              1.00000,        0.87500,        1.00000
+219,              1.00000,        0.87500,        0.87500
+220,              1.00000,        0.75000,        0.75000
+221,              1.00000,        0.75000,        0.87500
+222,              1.00000,        0.87500,        0.75000
+223,              1.00000,        0.87500,        0.62500
+224,              1.00000,        0.75000,        0.50000
+225,              1.00000,        0.75000,        0.62500
+226,              1.00000,        0.87500,        0.50000
+227,              1.00000,        0.50000,        1.00000
+228,              1.00000,        0.62500,        1.00000
+229,              1.00000,        0.62500,        0.87500
+230,              1.00000,        0.50000,        0.75000
+231,              1.00000,        0.50000,        0.87500
+232,              1.00000,        0.62500,        0.75000
+233,              1.00000,        0.62500,        0.62500
+234,              1.00000,        0.50000,        0.50000
+235,              1.00000,        0.50000,        0.62500
+236,              1.00000,        0.62500,        0.50000
+237,              1.00000,        0.25000,        1.00000
+238,              1.00000,        0.37500,        1.00000
+239,              1.00000,        0.37500,        0.87500
+240,              1.00000,        0.25000,        0.75000
+241,              1.00000,        0.25000,        0.87500
+242,              1.00000,        0.37500,        0.75000
+243,              1.00000,        0.37500,        0.62500
+244,              1.00000,        0.25000,        0.50000
+245,              1.00000,        0.25000,        0.62500
+246,              1.00000,        0.37500,        0.50000
+247,              1.00000,        0.12500,        1.00000
+248,              1.00000,        0.12500,        0.87500
+249,              1.00000,        0.12500,        0.75000
+250,              1.00000,        0.12500,        0.62500
+251,              1.00000,        0.12500,        0.50000
+252,              0.00000,        0.75000,        1.00000
+253,              0.00000,        0.87500,        1.00000
+254,              0.00000,        0.87500,        0.87500
+255,              0.00000,        0.75000,        0.75000
+256,              0.00000,        0.75000,        0.87500
+257,              0.00000,        0.87500,        0.75000
+258,              0.00000,        0.87500,        0.62500
+259,              0.00000,        0.75000,        0.50000
+260,              0.00000,        0.75000,        0.62500
+261,              0.00000,        0.87500,        0.50000
+262,              0.00000,        0.50000,        1.00000
+263,              0.00000,        0.62500,        1.00000
+264,              0.00000,        0.62500,        0.87500
+265,              0.00000,        0.50000,        0.75000
+266,              0.00000,        0.50000,        0.87500
+267,              0.00000,        0.62500,        0.75000
+268,              0.00000,        0.62500,        0.62500
+269,              0.00000,        0.50000,        0.50000
+270,              0.00000,        0.50000,        0.62500
+271,              0.00000,        0.62500,        0.50000
+272,              0.00000,        0.25000,        1.00000
+273,              0.00000,        0.37500,        1.00000
+274,              0.00000,        0.37500,        0.87500
+275,              0.00000,        0.25000,        0.75000
+276,              0.00000,        0.25000,        0.87500
+277,              0.00000,        0.37500,        0.75000
+278,              0.00000,        0.37500,        0.62500
+279,              0.00000,        0.25000,        0.50000
+280,              0.00000,        0.25000,        0.62500
+281,              0.00000,        0.37500,        0.50000
+282,              0.00000,        0.12500,        1.00000
+283,              0.00000,        0.12500,        0.87500
+284,              0.00000,        0.12500,        0.75000
+285,              0.00000,        0.12500,        0.62500
+286,              0.00000,        0.12500,        0.50000
+287,              0.87500,        0.87500,        1.00000
+288,              0.75000,        0.75000,        1.00000
+289,              0.87500,        0.75000,        1.00000
+290,              0.75000,        0.87500,        1.00000
+291,              0.62500,        0.87500,        1.00000
+292,              0.50000,        0.75000,        1.00000
+293,              0.62500,        0.75000,        1.00000
+294,              0.50000,        0.87500,        1.00000
+295,              0.37500,        0.87500,        1.00000
+296,              0.25000,        0.75000,        1.00000
+297,              0.37500,        0.75000,        1.00000
+298,              0.25000,        0.87500,        1.00000
+299,              0.12500,        0.87500,        1.00000
+300,              0.12500,        0.75000,        1.00000
+301,              0.87500,        0.62500,        1.00000
+302,              0.75000,        0.50000,        1.00000
+303,              0.87500,        0.50000,        1.00000
+304,              0.75000,        0.62500,        1.00000
+305,              0.62500,        0.62500,        1.00000
+306,              0.50000,        0.50000,        1.00000
+307,              0.62500,        0.50000,        1.00000
+308,              0.50000,        0.62500,        1.00000
+309,              0.37500,        0.62500,        1.00000
+310,              0.25000,        0.50000,        1.00000
+311,              0.37500,        0.50000,        1.00000
+312,              0.25000,        0.62500,        1.00000
+313,              0.12500,        0.62500,        1.00000
+314,              0.12500,        0.50000,        1.00000
+315,              0.87500,        0.37500,        1.00000
+316,              0.75000,        0.25000,        1.00000
+317,              0.87500,        0.25000,        1.00000
+318,              0.75000,        0.37500,        1.00000
+319,              0.62500,        0.37500,        1.00000
+320,              0.50000,        0.25000,        1.00000
+321,              0.62500,        0.25000,        1.00000
+322,              0.50000,        0.37500,        1.00000
+323,              0.37500,        0.37500,        1.00000
+324,              0.25000,        0.25000,        1.00000
+325,              0.37500,        0.25000,        1.00000
+326,              0.25000,        0.37500,        1.00000
+327,              0.12500,        0.37500,        1.00000
+328,              0.12500,        0.25000,        1.00000
+329,              0.87500,        0.12500,        1.00000
+330,              0.75000,        0.12500,        1.00000
+331,              0.62500,        0.12500,        1.00000
+332,              0.50000,        0.12500,        1.00000
+333,              0.37500,        0.12500,        1.00000
+334,              0.25000,        0.12500,        1.00000
+335,              0.12500,        0.12500,        1.00000
+336,              0.87500,        0.87500,        0.50000
+337,              0.75000,        0.75000,        0.50000
+338,              0.87500,        0.75000,        0.50000
+339,              0.75000,        0.87500,        0.50000
+340,              0.62500,        0.87500,        0.50000
+341,              0.50000,        0.75000,        0.50000
+342,              0.62500,        0.75000,        0.50000
+343,              0.50000,        0.87500,        0.50000
+344,              0.37500,        0.87500,        0.50000
+345,              0.25000,        0.75000,        0.50000
+346,              0.37500,        0.75000,        0.50000
+347,              0.25000,        0.87500,        0.50000
+348,              0.12500,        0.87500,        0.50000
+349,              0.12500,        0.75000,        0.50000
+350,              0.87500,        0.62500,        0.50000
+351,              0.75000,        0.50000,        0.50000
+352,              0.87500,        0.50000,        0.50000
+353,              0.75000,        0.62500,        0.50000
+354,              0.62500,        0.62500,        0.50000
+355,              0.50000,        0.50000,        0.50000
+356,              0.62500,        0.50000,        0.50000
+357,              0.50000,        0.62500,        0.50000
+358,              0.37500,        0.62500,        0.50000
+359,              0.25000,        0.50000,        0.50000
+360,              0.37500,        0.50000,        0.50000
+361,              0.25000,        0.62500,        0.50000
+362,              0.12500,        0.62500,        0.50000
+363,              0.12500,        0.50000,        0.50000
+364,              0.87500,        0.37500,        0.50000
+365,              0.75000,        0.25000,        0.50000
+366,              0.87500,        0.25000,        0.50000
+367,              0.75000,        0.37500,        0.50000
+368,              0.62500,        0.37500,        0.50000
+369,              0.50000,        0.25000,        0.50000
+370,              0.62500,        0.25000,        0.50000
+371,              0.50000,        0.37500,        0.50000
+372,              0.37500,        0.37500,        0.50000
+373,              0.25000,        0.25000,        0.50000
+374,              0.37500,        0.25000,        0.50000
+375,              0.25000,        0.37500,        0.50000
+376,              0.12500,        0.37500,        0.50000
+377,              0.12500,        0.25000,        0.50000
+378,              0.87500,        0.12500,        0.50000
+379,              0.75000,        0.12500,        0.50000
+380,              0.62500,        0.12500,        0.50000
+381,              0.50000,        0.12500,        0.50000
+382,              0.37500,        0.12500,        0.50000
+383,              0.25000,        0.12500,        0.50000
+384,              0.12500,        0.12500,        0.50000
+385,              0.83333,        0.33333,        0.12500
+386,              0.83333,        0.50000,        0.12500
+387,              0.83333,        0.83333,        0.25000
+388,              0.83333,        0.16667,        0.12500
+389,              0.66667,        0.16667,        0.12500
+390,              0.16667,        0.50000,        0.12500
+391,              0.16667,        0.66667,        0.12500
+392,              0.83333,        0.33333,        0.37500
+393,              0.83333,        0.50000,        0.37500
+394,              0.83333,        0.66667,        0.12500
+395,              0.33333,        0.16667,        0.12500
+396,              0.16667,        0.33333,        0.37500
+397,              0.16667,        0.50000,        0.37500
+398,              0.33333,        0.83333,        0.12500
+399,              0.50000,        0.83333,        0.12500
+400,              0.16667,        0.83333,        0.12500
+401,              0.83333,        0.16667,        0.37500
+402,              0.83333,        0.16667,        0.25000
+403,              0.50000,        0.83333,        0.37500
+404,              0.66667,        0.83333,        0.37500
+405,              0.16667,        0.33333,        0.12500
+406,              0.16667,        0.16667,        0.12500
+407,              0.16667,        0.16667,        0.25000
+408,              0.16667,        0.16667,        0.37500
+409,              0.83333,        0.66667,        0.37500
+410,              0.33333,        0.83333,        0.37500
+411,              0.50000,        0.50000,        0.25000
+412,              0.41667,        0.41667,        0.37500
+413,              0.25000,        0.41667,        0.25000
+414,              0.25000,        0.58333,        0.25000
+415,              0.33333,        0.16667,        0.37500
+416,              0.41667,        0.25000,        0.25000
+417,              0.50000,        0.16667,        0.37500
+418,              0.58333,        0.25000,        0.25000
+419,              0.66667,        0.16667,        0.37500
+420,              0.58333,        0.41667,        0.37500
+421,              0.41667,        0.58333,        0.37500
+422,              0.16667,        0.66667,        0.37500
+423,              0.16667,        0.83333,        0.25000
+424,              0.41667,        0.58333,        0.12500
+425,              0.41667,        0.75000,        0.25000
+426,              0.58333,        0.75000,        0.25000
+427,              0.58333,        0.58333,        0.37500
+428,              0.75000,        0.58333,        0.25000
+429,              0.75000,        0.41667,        0.25000
+430,              0.58333,        0.41667,        0.12500
+431,              0.50000,        0.16667,        0.12500
+432,              0.41667,        0.41667,        0.12500
+433,              0.58333,        0.58333,        0.12500
+434,              0.66667,        0.83333,        0.12500
+435,              0.16667,        0.83333,        0.37500
+436,              0.83333,        0.83333,        0.37500
+437,              0.83333,        0.83333,        0.12500
+438,              1.00000,        0.00000,        0.25000
+439,              1.00000,        0.00000,        0.50000
+440,              0.66667,        0.00000,        0.50000
+441,              1.00000,        0.00000,        0.37500
+442,              0.83333,        0.00000,        0.50000
+443,              0.83333,        0.00000,        0.37500
+444,              0.66667,        0.00000,        0.25000
+445,              0.66667,        0.00000,        0.37500
+446,              0.83333,        0.00000,        0.25000
+447,              1.00000,        0.00000,        0.00000
+448,              1.00000,        0.00000,        0.12500
+449,              0.83333,        0.00000,        0.12500
+450,              0.66667,        0.00000,        0.00000
+451,              0.66667,        0.00000,        0.12500
+452,              0.83333,        0.00000,        0.00000
+453,              0.33333,        0.00000,        0.50000
+454,              0.50000,        0.00000,        0.50000
+455,              0.50000,        0.00000,        0.37500
+456,              0.33333,        0.00000,        0.25000
+457,              0.33333,        0.00000,        0.37500
+458,              0.50000,        0.00000,        0.25000
+459,              0.50000,        0.00000,        0.12500
+460,              0.33333,        0.00000,        0.00000
+461,              0.33333,        0.00000,        0.12500
+462,              0.50000,        0.00000,        0.00000
+463,              0.00000,        0.00000,        0.50000
+464,              0.16667,        0.00000,        0.50000
+465,              0.16667,        0.00000,        0.37500
+466,              0.00000,        0.00000,        0.25000
+467,              0.00000,        0.00000,        0.37500
+468,              0.16667,        0.00000,        0.25000
+469,              0.16667,        0.00000,        0.12500
+470,              0.00000,        0.00000,        0.00000
+471,              0.00000,        0.00000,        0.12500
+472,              0.16667,        0.00000,        0.00000
+473,              0.66667,        1.00000,        0.50000
+474,              1.00000,        1.00000,        0.50000
+475,              1.00000,        1.00000,        0.25000
+476,              0.83333,        1.00000,        0.50000
+477,              1.00000,        1.00000,        0.37500
+478,              0.83333,        1.00000,        0.37500
+479,              0.66667,        1.00000,        0.25000
+480,              0.66667,        1.00000,        0.37500
+481,              0.83333,        1.00000,        0.25000
+482,              1.00000,        1.00000,        0.00000
+483,              1.00000,        1.00000,        0.12500
+484,              0.83333,        1.00000,        0.12500
+485,              0.66667,        1.00000,        0.00000
+486,              0.66667,        1.00000,        0.12500
+487,              0.83333,        1.00000,        0.00000
+488,              0.33333,        1.00000,        0.50000
+489,              0.50000,        1.00000,        0.50000
+490,              0.50000,        1.00000,        0.37500
+491,              0.33333,        1.00000,        0.25000
+492,              0.33333,        1.00000,        0.37500
+493,              0.50000,        1.00000,        0.25000
+494,              0.50000,        1.00000,        0.12500
+495,              0.33333,        1.00000,        0.00000
+496,              0.33333,        1.00000,        0.12500
+497,              0.50000,        1.00000,        0.00000
+498,              0.00000,        1.00000,        0.50000
+499,              0.16667,        1.00000,        0.50000
+500,              0.16667,        1.00000,        0.37500
+501,              0.00000,        1.00000,        0.25000
+502,              0.00000,        1.00000,        0.37500
+503,              0.16667,        1.00000,        0.25000
+504,              0.16667,        1.00000,        0.12500
+505,              0.00000,        1.00000,        0.00000
+506,              0.00000,        1.00000,        0.12500
+507,              0.16667,        1.00000,        0.00000
+508,              1.00000,        0.66667,        0.50000
+509,              1.00000,        0.83333,        0.50000
+510,              1.00000,        0.83333,        0.37500
+511,              1.00000,        0.66667,        0.25000
+512,              1.00000,        0.66667,        0.37500
+513,              1.00000,        0.83333,        0.25000
+514,              1.00000,        0.83333,        0.12500
+515,              1.00000,        0.66667,        0.00000
+516,              1.00000,        0.66667,        0.12500
+517,              1.00000,        0.83333,        0.00000
+518,              1.00000,        0.33333,        0.50000
+519,              1.00000,        0.50000,        0.50000
+520,              1.00000,        0.50000,        0.37500
+521,              1.00000,        0.33333,        0.25000
+522,              1.00000,        0.33333,        0.37500
+523,              1.00000,        0.50000,        0.25000
+524,              1.00000,        0.50000,        0.12500
+525,              1.00000,        0.33333,        0.00000
+526,              1.00000,        0.33333,        0.12500
+527,              1.00000,        0.50000,        0.00000
+528,              1.00000,        0.16667,        0.50000
+529,              1.00000,        0.16667,        0.37500
+530,              1.00000,        0.16667,        0.25000
+531,              1.00000,        0.16667,        0.12500
+532,              1.00000,        0.16667,        0.00000
+533,              0.00000,        0.66667,        0.50000
+534,              0.00000,        0.83333,        0.50000
+535,              0.00000,        0.83333,        0.37500
+536,              0.00000,        0.66667,        0.25000
+537,              0.00000,        0.66667,        0.37500
+538,              0.00000,        0.83333,        0.25000
+539,              0.00000,        0.83333,        0.12500
+540,              0.00000,        0.66667,        0.00000
+541,              0.00000,        0.66667,        0.12500
+542,              0.00000,        0.83333,        0.00000
+543,              0.00000,        0.33333,        0.50000
+544,              0.00000,        0.50000,        0.50000
+545,              0.00000,        0.50000,        0.37500
+546,              0.00000,        0.33333,        0.25000
+547,              0.00000,        0.33333,        0.37500
+548,              0.00000,        0.50000,        0.25000
+549,              0.00000,        0.50000,        0.12500
+550,              0.00000,        0.33333,        0.00000
+551,              0.00000,        0.33333,        0.12500
+552,              0.00000,        0.50000,        0.00000
+553,              0.00000,        0.16667,        0.50000
+554,              0.00000,        0.16667,        0.37500
+555,              0.00000,        0.16667,        0.25000
+556,              0.00000,        0.16667,        0.12500
+557,              0.00000,        0.16667,        0.00000
+558,              0.83333,        0.83333,        0.50000
+559,              0.66667,        0.66667,        0.50000
+560,              0.83333,        0.66667,        0.50000
+561,              0.66667,        0.83333,        0.50000
+562,              0.50000,        0.83333,        0.50000
+563,              0.33333,        0.66667,        0.50000
+564,              0.50000,        0.66667,        0.50000
+565,              0.33333,        0.83333,        0.50000
+566,              0.16667,        0.83333,        0.50000
+567,              0.16667,        0.66667,        0.50000
+568,              0.83333,        0.50000,        0.50000
+569,              0.66667,        0.33333,        0.50000
+570,              0.83333,        0.33333,        0.50000
+571,              0.66667,        0.50000,        0.50000
+572,              0.50000,        0.50000,        0.50000
+573,              0.33333,        0.33333,        0.50000
+574,              0.50000,        0.33333,        0.50000
+575,              0.33333,        0.50000,        0.50000
+576,              0.16667,        0.50000,        0.50000
+577,              0.16667,        0.33333,        0.50000
+578,              0.83333,        0.16667,        0.50000
+579,              0.66667,        0.16667,        0.50000
+580,              0.50000,        0.16667,        0.50000
+581,              0.33333,        0.16667,        0.50000
+582,              0.16667,        0.16667,        0.50000
+583,              0.83333,        0.83333,        0.00000
+584,              0.66667,        0.66667,        0.00000
+585,              0.83333,        0.66667,        0.00000
+586,              0.66667,        0.83333,        0.00000
+587,              0.50000,        0.83333,        0.00000
+588,              0.33333,        0.66667,        0.00000
+589,              0.50000,        0.66667,        0.00000
+590,              0.33333,        0.83333,        0.00000
+591,              0.16667,        0.83333,        0.00000
+592,              0.16667,        0.66667,        0.00000
+593,              0.83333,        0.50000,        0.00000
+594,              0.66667,        0.33333,        0.00000
+595,              0.83333,        0.33333,        0.00000
+596,              0.66667,        0.50000,        0.00000
+597,              0.50000,        0.50000,        0.00000
+598,              0.33333,        0.33333,        0.00000
+599,              0.50000,        0.33333,        0.00000
+600,              0.33333,        0.50000,        0.00000
+601,              0.16667,        0.50000,        0.00000
+602,              0.16667,        0.33333,        0.00000
+603,              0.83333,        0.16667,        0.00000
+604,              0.66667,        0.16667,        0.00000
+605,              0.50000,        0.16667,        0.00000
+606,              0.33333,        0.16667,        0.00000
+607,              0.16667,        0.16667,        0.00000
+**
+**ELSET COUNT = 177
+**HWCOLOR COMP         54     0
+*ELEMENT, TYPE=C3D10, ELSET=UPPER
+         1,       292,       187,       288,         1,       294,       291,       293,
+         2,         3,         4
+         2,       302,       227,       230,         1,       303,       231,         5,
+         6,         7,         8
+         3,       302,       288,       227,         1,       304,       301,       303,
+         6,         4,         7
+         4,       227,       288,       217,         1,       301,       289,       228,
+         7,         4,         9
+         5,       227,       217,       220,         1,       228,       221,       229,
+         7,         9,        10
+         6,       227,       220,       230,         1,       229,       232,       231,
+         7,        10,         8
+         7,       337,       230,       220,         1,        11,       232,        12,
+        13,         8,        10
+         8,       351,       230,       337,        14,        15,        11,       353,
+        16,        17,        18
+         9,       341,       337,       190,         1,       342,        19,        20,
+        21,        13,        22
+        10,       337,       178,       190,         1,        23,       192,        19,
+        13,        24,        22
+        11,       337,       220,       178,         1,        12,        25,        23,
+        13,        10,        24
+        12,       217,       178,       220,         1,        26,        25,       221,
+         9,        24,        10
+        13,       217,       172,       178,         1,       287,       179,        26,
+         9,        27,        24
+        14,       288,       172,       217,         1,       290,       287,       289,
+         4,        27,         9
+        15,       288,       187,       172,         1,       291,       188,       290,
+         4,         3,        27
+        16,       187,       178,       172,         1,       189,       179,       188,
+         3,        24,        27
+        17,       306,       292,       302,        14,       308,       305,       307,
+        28,        29,        30
+        18,       320,       306,       316,        14,       322,       319,       321,
+        31,        28,        32
+        19,       316,       306,       302,        14,       319,       307,       318,
+        32,        28,        30
+        20,       316,       302,       240,        14,       318,        33,        34,
+        32,        30,        35
+        21,       369,       365,       355,        14,       370,       368,       371,
+        36,        37,        38
+        22,       365,       351,       355,        14,       367,       356,       368,
+        37,        16,        38
+        23,       365,       240,       351,        14,        39,        40,       367,
+        37,        35,        16
+        24,       351,       240,       230,        14,        40,       242,        15,
+        16,        35,        17
+        25,       302,       230,       240,        14,         5,       242,        33,
+        30,        17,        35
+        26,       302,         1,       230,        14,         6,         8,         5,
+        30,        41,        17
+        27,       302,       292,       288,         1,       305,       293,       304,
+         6,         2,         4
+        28,       302,       292,         1,        14,       305,         2,         6,
+        30,        29,        41
+        29,       337,       224,       184,       178,       338,       336,       339,
+        23,        42,       185
+        30,       355,       351,       341,        14,       356,       354,       357,
+        38,        16,        43
+        31,       351,       337,       341,        14,       353,       342,       354,
+        16,        18,        43
+        32,       337,       230,         1,        14,        11,         8,        13,
+        18,        17,        41
+        33,       365,       136,       244,       133,       378,       251,       366,
+        44,       138,        45
+        34,       341,       337,         1,        14,       342,        13,        21,
+        43,        18,        41
+        35,       129,       240,       127,       133,        46,       249,       132,
+       134,        47,       135
+        36,       292,         1,        14,        48,         2,        41,        29,
+        49,        50,        51
+        37,       306,       292,        14,        48,       308,        29,        28,
+        52,        49,        51
+        38,       306,       296,       292,        48,       309,       297,       308,
+        52,        53,        49
+        39,       310,       296,       306,        48,       312,       309,       311,
+        54,        53,        52
+        40,       310,       255,       296,        48,        55,        56,       312,
+        54,        57,        53
+        41,       310,       265,       255,        48,        58,       267,        55,
+        54,        59,        57
+        42,       355,        14,       341,        48,        38,        43,       357,
+        60,        51,        61
+        43,       341,        14,         1,        48,        43,        41,        21,
+        61,        51,        50
+        44,       341,         1,       190,        48,        21,        22,        20,
+        61,        50,        62
+        45,       341,       190,       194,        48,        20,       195,       343,
+        61,        62,        63
+        46,       341,       194,       204,        48,       343,       206,       344,
+        61,        63,        64
+        47,       345,       341,       204,        48,       346,       344,       347,
+        65,        61,        64
+        48,       355,       341,       345,        48,       357,       346,       358,
+        60,        61,        65
+        49,       359,       355,       345,        48,       360,       358,       361,
+        66,        60,        65
+        50,       359,       345,       265,        48,       361,        67,        68,
+        66,        65,        59
+        51,       345,       255,       265,        48,        69,       267,        67,
+        65,        57,        59
+        52,       259,       345,       214,       255,       349,       348,       261,
+       260,        69,       258
+        53,       345,       204,       255,        48,       347,        70,        69,
+        65,        64,        57
+        54,       204,       194,       200,        48,       206,       203,       205,
+        64,        63,        71
+        55,       194,       190,       200,        48,       195,       202,       203,
+        63,        62,        71
+        56,       197,       200,       190,        48,       201,       202,       199,
+        72,        71,        62
+        57,       296,       200,       197,        48,        73,       201,       298,
+        53,        71,        72
+        58,       296,       197,       292,        48,       298,       295,       297,
+        53,        72,        49
+        59,       292,       197,       187,        48,       295,       198,       294,
+        49,        72,        74
+        60,       292,       187,         1,        48,       294,         3,         2,
+        49,        74,        50
+        61,       197,       190,       187,        48,       199,       191,       198,
+        72,        62,        74
+        62,       296,       255,       200,        48,        56,        75,        73,
+        53,        57,        71
+        63,       142,       145,       320,        76,       146,        77,       332,
+        78,        79,        80
+        64,       142,       133,       145,        76,       144,       147,       146,
+        78,        81,        79
+        65,       142,       129,       133,        76,       143,       134,       144,
+        78,        82,        81
+        66,       142,       320,       129,        76,       332,       331,       143,
+        78,        80,        82
+        67,       320,       316,       129,        76,       321,       330,       331,
+        80,        83,        82
+        68,       320,        14,       316,        76,        31,        32,       321,
+        80,        84,        83
+        69,       316,        14,       240,        76,        32,        35,        34,
+        83,        84,        85
+        70,       129,       240,       133,        76,        46,        47,       134,
+        82,        85,        81
+        71,       129,       316,       240,        76,       330,        34,        46,
+        82,        83,        85
+        72,       365,       133,       240,        76,        44,        47,        39,
+        86,        81,        85
+        73,       365,       240,        14,        76,        39,        35,        37,
+        86,        85,        84
+        74,       369,       365,        14,        76,       370,        37,        36,
+        87,        86,        84
+        75,       369,       139,       365,        76,       380,       379,       370,
+        87,        88,        86
+        76,       149,       139,       369,        76,       151,       380,       381,
+        89,        88,        87
+        77,       149,       145,       139,        76,       150,       148,       151,
+        89,        79,        88
+        78,       373,       149,       369,        90,       382,       381,       374,
+        91,        92,        93
+        79,       139,       145,       133,        76,       148,       147,       140,
+        88,        79,        81
+        80,       139,       133,       365,        76,       140,        44,       379,
+        88,        81,        86
+        81,       310,        48,       306,        94,        54,        52,       311,
+        95,        96,        97
+        82,       310,       265,        48,        94,        58,        59,        54,
+        95,        98,        96
+        83,       324,       310,       320,        94,       326,       323,       325,
+        99,        95,       100
+        84,       373,       369,       359,        94,       374,       372,       375,
+       101,       102,       103
+        85,       369,        76,        14,        94,        87,        84,        36,
+       102,       104,       105
+        86,       320,        14,        76,        94,        31,        84,        80,
+       100,       105,       104
+        87,       320,       306,        14,        94,       322,        28,        31,
+       100,        97,       105
+        88,       320,       310,       306,        94,       323,       311,       322,
+       100,        95,        97
+        89,       306,        48,        14,        94,        52,        51,        28,
+        97,        96,       105
+        90,       355,        14,        48,        94,        38,        51,        60,
+       106,       105,        96
+        91,       369,        14,       355,        94,        36,        38,       371,
+       102,       105,       106
+        92,       369,       355,       359,        94,       371,       360,       372,
+       102,       106,       103
+        93,       359,       355,        48,        94,       360,        60,        66,
+       103,       106,        96
+        94,       359,        48,       265,        94,        66,        59,        68,
+       103,        96,        98
+        95,       324,       145,       155,        90,       107,       157,       108,
+       109,       110,       111
+        96,       324,       155,       275,        90,       108,       112,       113,
+       109,       111,       114
+        97,       324,       275,       265,        90,       113,       277,       115,
+       109,       114,       116
+        98,       324,       265,       310,        94,       115,        58,       326,
+        99,        98,        95
+        99,       224,       181,       184,       178,       226,       186,       336,
+        42,       183,       185
+       100,       324,        94,       320,        90,        99,       100,       325,
+       109,       117,       118
+       101,       324,       320,       145,        90,       325,        77,       107,
+       109,       118,       110
+       102,       320,        76,       145,        90,        80,        79,        77,
+       118,       119,       110
+       103,       320,        94,        76,        90,       100,       104,        80,
+       118,       117,       119
+       104,       373,       369,        94,        90,       374,       102,       101,
+        91,        93,       117
+       105,       149,       145,        76,        90,       150,        79,        89,
+        92,       110,       119
+       106,       162,       275,       152,       155,       283,       120,       163,
+       164,       112,       156
+       107,       149,        76,       369,        90,        89,        87,       381,
+        92,       119,        93
+       108,       159,       149,       373,        90,       161,       382,       383,
+       121,        92,        91
+       109,       159,       155,       149,        90,       160,       158,       161,
+       121,       111,        92
+       110,       159,       275,       155,        90,       122,       112,       160,
+       121,       114,       111
+       111,       279,       275,       159,        90,       280,       122,       384,
+       123,       114,       121
+       112,       279,       159,       373,        90,       384,       383,       377,
+       123,       121,        91
+       113,       279,       373,       269,        90,       377,       376,       281,
+       123,        91,       124
+       114,       279,       269,       275,        90,       281,       278,       280,
+       123,       124,       114
+       115,       269,       265,       275,        90,       270,       277,       278,
+       124,       116,       114
+       116,       269,       359,       265,        90,       363,        68,       270,
+       124,       125,       116
+       117,       373,       359,       269,        90,       375,       363,       376,
+        91,       125,       124
+       118,       373,        94,       359,        90,       101,       103,       375,
+        91,       117,       125
+       119,       359,        94,       265,        90,       103,        98,        68,
+       125,       117,       116
+       120,       324,       265,        94,        90,       115,        98,        99,
+       109,       116,       117
+       121,       369,        76,        94,        90,        87,       104,       102,
+        93,       119,       117
+       122,       149,       155,       145,        90,       158,       157,       150,
+        92,       111,       110
+       123,       217,       174,       220,       178,       219,       222,       221,
+        26,       180,        25
+       124,       345,       204,       214,       255,       347,       216,       348,
+        69,        70,       258
+       125,       337,       184,       190,       178,       339,       193,        19,
+        23,       185,       192
+       126,       337,       220,       224,       178,        12,       225,       338,
+        23,        25,        42
+       127,       217,       172,       174,       178,       287,       177,       219,
+        26,       179,       180
+       128,       204,       210,       255,       200,       213,       257,        70,
+       205,       212,        75
+       129,       234,       224,       337,       230,       236,       338,       350,
+       235,       233,        11
+       130,       181,       220,       174,       178,       223,       222,       182,
+       183,        25,       180
+       131,       337,       230,       224,       220,        11,       233,       338,
+        12,       232,       225
+       132,       341,       337,       194,       190,       342,       340,       343,
+        20,        19,       195
+       133,       237,       240,       302,       230,       241,        33,       315,
+       239,       242,         5
+       134,       359,       345,       259,       265,       361,       349,       362,
+        68,        67,       268
+       135,       237,       302,       227,       230,       315,       303,       238,
+       239,         5,       231
+       136,       351,       240,       234,       230,        40,       243,       352,
+        15,       242,       235
+       137,       351,       234,       337,       230,       352,       350,       353,
+        15,       235,        11
+       138,       244,       127,       240,       133,       250,       249,       245,
+        45,       135,        47
+       139,       136,       127,       244,       133,       137,       250,       251,
+       138,       135,        45
+       140,       337,       184,       194,       190,       339,       196,       340,
+        19,       193,       195
+       141,       244,       234,       351,       240,       246,       352,       364,
+       245,       243,        40
+       142,       365,       244,       351,       240,       366,       364,       367,
+        39,       245,        40
+       143,       152,       155,       324,       145,       156,       108,       334,
+       154,       157,       107
+       144,       316,       237,       128,       240,       317,       247,       329,
+        34,       241,       248
+       145,       129,       316,       128,       240,       330,       329,       131,
+        46,        34,       248
+       146,       252,       210,       207,       200,       254,       211,       253,
+       126,       212,       209
+       147,       269,       359,       259,       265,       363,       362,       271,
+       270,        68,       268
+       148,       324,       262,       310,       265,       327,       314,       326,
+       115,       266,        58
+       149,       310,       252,       296,       255,       313,       300,       312,
+        55,       256,        56
+       150,       324,       320,       142,       145,       325,       332,       333,
+       107,        77,       146
+       151,       162,       165,       275,       155,       166,       284,       283,
+       164,       167,       112
+       152,       279,       169,       159,       165,       286,       171,       384,
+       285,       170,       168
+       153,       214,       255,       204,       210,       258,        70,       216,
+       215,       257,       213
+       154,       217,       172,       173,       174,       287,       175,       218,
+       219,       177,       176
+       155,       272,       262,       324,       265,       273,       327,       328,
+       274,       266,       115
+       156,       316,       302,       237,       240,       318,       315,       317,
+        34,        33,       241
+       157,       129,       128,       127,       240,       131,       130,       132,
+        46,       248,       249
+       158,       139,       136,       365,       133,       141,       378,       379,
+       140,       138,        44
+       159,       159,       275,       165,       155,       122,       284,       168,
+       160,       112,       167
+       160,       152,       272,       324,       275,       335,       328,       334,
+       120,       276,       113
+       161,       259,       265,       345,       255,       268,        67,       349,
+       260,       267,        69
+       162,       262,       252,       310,       255,       263,       313,       314,
+       264,       256,        55
+       163,       252,       296,       255,       200,       300,        56,       256,
+       126,        73,        75
+       164,       204,       200,       255,        48,       205,        75,        70,
+        64,        71,        57
+       165,       365,       244,       240,       133,       366,       245,        39,
+        44,        45,        47
+       166,       279,       165,       159,       275,       285,       168,       384,
+       280,       284,       122
+       167,       224,       220,       181,       178,       225,       223,       226,
+        42,        25,       183
+       168,       296,       207,       197,       200,       299,       208,       298,
+        73,       209,       201
+       169,       272,       324,       275,       265,       328,       113,       276,
+       274,       115,       277
+       170,       252,       207,       296,       200,       253,       299,       300,
+       126,       209,        73
+       171,       252,       255,       210,       200,       256,       257,       254,
+       126,        75,       212
+       172,       162,       272,       152,       275,       282,       335,       163,
+       283,       276,       120
+       173,       262,       310,       265,       255,       314,        58,       266,
+       264,        55,       267
+       174,       152,       275,       324,       155,       120,       113,       334,
+       156,       112,       108
+       175,       152,       324,       142,       145,       334,       333,       153,
+       154,       107,       146
+       176,       187,       190,       178,         1,       191,       192,       189,
+         3,        22,        24
+       177,       187,       190,         1,        48,       191,        22,         3,
+        74,        62,        50
+**
+**ELSET COUNT = 92
+**HWCOLOR COMP          1     0
+*ELEMENT, TYPE=C3D10, ELSET=LOWER
+       354,       594,       525,       584,       521,       595,       593,       596,
+       385,       526,       386
+       355,       482,       511,       475,       479,       514,       513,       483,
+       484,       387,       481
+       356,       594,       447,       438,       444,       603,       448,       388,
+       389,       449,       446
+       357,       540,       546,       588,       536,       549,       390,       592,
+       541,       548,       391
+       358,       569,       518,       521,       511,       570,       522,       392,
+       393,       520,       523
+       359,       584,       521,       515,       511,       386,       524,       585,
+       394,       523,       516
+       360,       460,       450,       598,       456,       462,       605,       606,
+       461,       459,       395
+       361,       543,       573,       546,       536,       577,       396,       547,
+       545,       397,       548
+       362,       588,       584,       495,       491,       589,       587,       590,
+       398,       399,       496
+       363,       588,       495,       501,       491,       590,       504,       400,
+       398,       496,       503
+       364,       588,       495,       505,       501,       590,       507,       591,
+       400,       504,       506
+       365,       588,       505,       536,       501,       591,       539,       391,
+       400,       506,       538
+       366,       450,       447,       594,       444,       452,       603,       604,
+       451,       449,       389
+       367,       440,       521,       438,       444,       401,       530,       443,
+       445,       402,       446
+       368,       563,       488,       559,       479,       565,       562,       564,
+       403,       490,       404
+       369,       550,       598,       546,       456,       602,       405,       551,
+       406,       395,       407
+       370,       463,       546,       543,       456,       554,       547,       553,
+       465,       407,       408
+       371,       440,       569,       439,       521,       579,       578,       442,
+       401,       392,       529
+       372,       569,       559,       518,       511,       571,       568,       570,
+       393,       409,       520
+       373,       463,       543,       453,       456,       553,       582,       464,
+       465,       408,       457
+       374,       550,       466,       460,       456,       556,       469,       607,
+       406,       468,       461
+       375,       563,       491,       488,       479,       410,       492,       565,
+       403,       493,       490
+       376,       573,       546,       536,       411,       396,       548,       397,
+       412,       413,       414
+       377,       573,       456,       546,       411,       415,       407,       396,
+       412,       416,       413
+       378,       573,       444,       456,       411,       417,       458,       415,
+       412,       418,       416
+       379,       573,       569,       444,       411,       574,       419,       417,
+       412,       420,       418
+       380,       573,       563,       569,       411,       575,       572,       574,
+       412,       421,       420
+       381,       573,       536,       563,       411,       397,       422,       575,
+       412,       414,       421
+       382,       588,       491,       536,       411,       398,       423,       391,
+       424,       425,       414
+       383,       563,       536,       491,       411,       422,       423,       410,
+       421,       414,       425
+       384,       563,       491,       479,       411,       410,       493,       403,
+       421,       425,       426
+       385,       563,       479,       559,       411,       403,       404,       564,
+       421,       426,       427
+       386,       569,       563,       559,       411,       572,       564,       571,
+       420,       421,       427
+       387,       569,       559,       511,       411,       571,       409,       393,
+       420,       427,       428
+       388,       569,       511,       521,       411,       393,       523,       392,
+       420,       428,       429
+       389,       594,       444,       521,       411,       389,       402,       385,
+       430,       418,       429
+       390,       569,       521,       444,       411,       392,       402,       419,
+       420,       429,       418
+       391,       594,       456,       444,       411,       431,       458,       389,
+       430,       416,       418
+       392,       598,       456,       594,       411,       395,       431,       599,
+       432,       416,       430
+       393,       598,       594,       588,       411,       599,       597,       600,
+       432,       430,       424
+       394,       598,       588,       546,       411,       600,       390,       405,
+       432,       424,       413
+       395,       598,       546,       456,       411,       405,       407,       395,
+       432,       413,       416
+       396,       588,       536,       546,       411,       391,       548,       390,
+       424,       414,       413
+       397,       594,       584,       588,       411,       596,       589,       597,
+       430,       433,       424
+       398,       594,       521,       584,       411,       385,       386,       596,
+       430,       429,       433
+       399,       584,       521,       511,       411,       386,       523,       394,
+       433,       429,       428
+       400,       559,       479,       511,       411,       404,       387,       409,
+       427,       426,       428
+       401,       584,       511,       479,       411,       394,       387,       434,
+       433,       428,       426
+       402,       584,       479,       491,       411,       434,       493,       399,
+       433,       426,       425
+       403,       588,       584,       491,       411,       589,       399,       398,
+       424,       433,       425
+       404,       533,       498,       563,       491,       534,       566,       567,
+       435,       500,       410
+       405,       559,       508,       511,       479,       560,       512,       409,
+       404,       436,       387
+       406,       515,       511,       482,       479,       516,       514,       517,
+       437,       387,       484
+       407,       594,       438,       525,       521,       388,       531,       595,
+       385,       530,       526
+       408,       594,       447,       525,       438,       603,       532,       595,
+       388,       448,       531
+       409,       598,       450,       594,       456,       605,       604,       599,
+       395,       459,       431
+       410,       508,       473,       475,       479,       558,       478,       510,
+       436,       480,       481
+       411,       515,       482,       485,       479,       517,       487,       583,
+       437,       484,       486
+       412,       508,       473,       474,       475,       558,       476,       509,
+       510,       478,       477
+       413,       559,       473,       508,       479,       561,       558,       560,
+       404,       480,       436
+       414,       518,       559,       508,       511,       568,       560,       519,
+       520,       409,       512
+       415,       563,       498,       488,       491,       566,       499,       565,
+       410,       500,       492
+       416,       559,       488,       473,       479,       562,       489,       561,
+       404,       490,       480
+       417,       584,       515,       485,       479,       585,       583,       586,
+       434,       437,       486
+       418,       588,       501,       536,       491,       400,       538,       391,
+       398,       503,       423
+       419,       594,       438,       521,       444,       388,       530,       385,
+       389,       446,       402
+       420,       533,       501,       498,       491,       535,       502,       534,
+       435,       503,       500
+       421,       543,       533,       573,       536,       544,       576,       577,
+       545,       537,       397
+       422,       569,       518,       439,       521,       570,       528,       578,
+       392,       522,       529
+       423,       533,       563,       536,       491,       567,       422,       537,
+       435,       410,       423
+       424,       440,       569,       521,       444,       579,       392,       401,
+       445,       419,       402
+       425,       450,       594,       456,       444,       604,       431,       459,
+       451,       389,       458
+       426,       573,       569,       440,       444,       574,       579,       580,
+       417,       419,       445
+       427,       550,       546,       466,       456,       551,       555,       556,
+       406,       407,       468
+       428,       440,       439,       438,       521,       442,       441,       443,
+       401,       529,       530
+       429,       508,       475,       511,       479,       510,       513,       512,
+       436,       481,       387
+       430,       453,       573,       440,       444,       581,       580,       454,
+       455,       417,       445
+       431,       453,       456,       573,       444,       457,       415,       581,
+       455,       458,       417
+       432,       550,       460,       598,       456,       607,       606,       602,
+       406,       461,       395
+       433,       453,       543,       573,       456,       582,       577,       581,
+       457,       408,       415
+       434,       543,       546,       573,       456,       547,       396,       577,
+       408,       407,       415
+       435,       540,       588,       505,       536,       592,       591,       542,
+       541,       391,       539
+       436,       573,       533,       563,       536,       576,       567,       575,
+       397,       537,       422
+       437,       550,       598,       540,       546,       602,       601,       552,
+       551,       405,       549
+       438,       533,       536,       501,       491,       537,       538,       535,
+       435,       423,       503
+       439,       550,       470,       460,       466,       557,       472,       607,
+       556,       471,       469
+       440,       598,       588,       540,       546,       600,       592,       601,
+       405,       390,       549
+       441,       584,       485,       491,       479,       586,       494,       399,
+       434,       486,       493
+       442,       584,       485,       495,       491,       586,       497,       587,
+       399,       494,       496
+       443,       584,       511,       515,       479,       394,       516,       585,
+       434,       387,       437
+       444,       525,       515,       584,       521,       527,       585,       593,
+       526,       524,       386
+       445,       463,       466,       546,       456,       467,       555,       554,
+       465,       468,       407
+**
+**ELSET COUNT = 18
+*SURFACE,TYPE=ELEMENT,NAME=LOWER_TO_UPPER
+373,S1
+433,S1
+430,S1
+426,S1
+371,S1
+422,S1
+421,S1
+436,S1
+380,S1
+386,S1
+372,S1
+414,S1
+404,S1
+415,S1
+368,S1
+416,S1
+413,S1
+412,S1
+**
+**ELSET COUNT = 18
+*SURFACE,TYPE=ELEMENT,NAME=LOWER_BOTTOM
+439,S1
+432,S1
+360,S1
+409,S1
+366,S1
+408,S1
+437,S1
+440,S1
+393,S1
+397,S1
+354,S1
+444,S1
+435,S1
+364,S1
+362,S1
+442,S1
+417,S1
+411,S1
+**
+**ELSET COUNT = 12
+*SURFACE,TYPE=ELEMENT,NAME=LOWER_SYM13
+428,S1
+445,S2
+360,S2
+431,S2
+366,S2
+439,S3
+374,S3
+356,S3
+373,S4
+425,S4
+430,S4
+367,S4
+**
+**ELSET COUNT = 12
+*SURFACE,TYPE=ELEMENT,NAME=LOWER_SYM23
+427,S1
+445,S1
+370,S1
+438,S1
+420,S1
+439,S2
+357,S2
+421,S2
+365,S3
+437,S4
+361,S4
+435,S4
+**
+**ELSET COUNT = 32
+*SURFACE,TYPE=ELEMENT,NAME=UPPER_TOP
+172,S1
+160,S1
+175,S1
+150,S1
+66,S1
+67,S1
+145,S1
+144,S1
+155,S1
+148,S1
+83,S1
+88,S1
+18,S1
+19,S1
+156,S1
+135,S1
+162,S1
+149,S1
+39,S1
+38,S1
+17,S1
+27,S1
+3,S1
+4,S1
+170,S1
+168,S1
+58,S1
+59,S1
+1,S1
+15,S1
+14,S1
+154,S1
+**
+**ELSET COUNT = 32
+*SURFACE,TYPE=ELEMENT,NAME=UPPER_TO_LOWER
+152,S1
+112,S1
+108,S1
+78,S1
+76,S1
+75,S1
+158,S1
+33,S1
+113,S1
+117,S1
+84,S1
+92,S1
+21,S1
+22,S1
+142,S1
+141,S1
+147,S1
+134,S1
+49,S1
+48,S1
+30,S1
+31,S1
+137,S1
+129,S1
+52,S1
+124,S1
+47,S1
+46,S1
+132,S1
+140,S1
+29,S1
+99,S1
+**
+**ELSET COUNT = 16
+*SURFACE,TYPE=ELEMENT,NAME=UPPER_SYM13
+109,S1
+122,S1
+77,S1
+79,S1
+64,S1
+65,S1
+157,S1
+151,S2
+143,S2
+158,S2
+139,S2
+152,S3
+159,S4
+106,S4
+175,S4
+35,S4
+**
+**ELSET COUNT = 16
+*SURFACE,TYPE=ELEMENT,NAME=UPPER_SYM23
+151,S1
+114,S1
+115,S1
+171,S1
+146,S1
+152,S2
+166,S2
+172,S2
+155,S2
+161,S2
+162,S2
+153,S2
+169,S4
+147,S4
+173,S4
+52,S4
+**
+**Property Definitions
+**
+*SOLID SECTION, ELSET=UPPER, MATERIAL=Def_Material
+*SOLID SECTION, ELSET=LOWER, MATERIAL=Def_Material
+**
+**Material Definitions
+**
+**Material:Def_Material
+*MATERIAL,NAME=Def_Material
+*ELASTIC,TYPE=ISO
+2.08000e+005,3.00000e-001
+*DENSITY
+7.80000e-009,
+*SPECIFIC HEAT
+5.00000e-001
+*CONDUCTIVITY
+4.98100e-002
+**
\ No newline at end of file
diff --git a/test/test_problems_contact_3d/tet4.inp b/test/test_problems_contact_3d/tet4.inp
new file mode 100644
index 0000000..fabfbb8
--- /dev/null
+++ b/test/test_problems_contact_3d/tet4.inp
@@ -0,0 +1,589 @@
+**NSET COUNT = 116
+*NODE
+127,              1.00000,        0.00000,        0.75000
+128,              1.00000,        0.00000,        1.00000
+129,              0.75000,        0.00000,        1.00000
+133,              0.75000,        0.00000,        0.75000
+136,              1.00000,        0.00000,        0.50000
+139,              0.75000,        0.00000,        0.50000
+142,              0.50000,        0.00000,        1.00000
+145,              0.50000,        0.00000,        0.75000
+149,              0.50000,        0.00000,        0.50000
+152,              0.25000,        0.00000,        1.00000
+155,              0.25000,        0.00000,        0.75000
+159,              0.25000,        0.00000,        0.50000
+162,              0.00000,        0.00000,        1.00000
+165,              0.00000,        0.00000,        0.75000
+169,              0.00000,        0.00000,        0.50000
+172,              0.75000,        1.00000,        1.00000
+173,              1.00000,        1.00000,        1.00000
+174,              1.00000,        1.00000,        0.75000
+178,              0.75000,        1.00000,        0.75000
+181,              1.00000,        1.00000,        0.50000
+184,              0.75000,        1.00000,        0.50000
+187,              0.50000,        1.00000,        1.00000
+190,              0.50000,        1.00000,        0.75000
+194,              0.50000,        1.00000,        0.50000
+197,              0.25000,        1.00000,        1.00000
+200,              0.25000,        1.00000,        0.75000
+204,              0.25000,        1.00000,        0.50000
+207,              0.00000,        1.00000,        1.00000
+210,              0.00000,        1.00000,        0.75000
+214,              0.00000,        1.00000,        0.50000
+217,              1.00000,        0.75000,        1.00000
+220,              1.00000,        0.75000,        0.75000
+224,              1.00000,        0.75000,        0.50000
+227,              1.00000,        0.50000,        1.00000
+230,              1.00000,        0.50000,        0.75000
+234,              1.00000,        0.50000,        0.50000
+237,              1.00000,        0.25000,        1.00000
+240,              1.00000,        0.25000,        0.75000
+244,              1.00000,        0.25000,        0.50000
+252,              0.00000,        0.75000,        1.00000
+255,              0.00000,        0.75000,        0.75000
+259,              0.00000,        0.75000,        0.50000
+262,              0.00000,        0.50000,        1.00000
+265,              0.00000,        0.50000,        0.75000
+269,              0.00000,        0.50000,        0.50000
+272,              0.00000,        0.25000,        1.00000
+275,              0.00000,        0.25000,        0.75000
+279,              0.00000,        0.25000,        0.50000
+288,              0.75000,        0.75000,        1.00000
+292,              0.50000,        0.75000,        1.00000
+296,              0.25000,        0.75000,        1.00000
+302,              0.75000,        0.50000,        1.00000
+306,              0.50000,        0.50000,        1.00000
+310,              0.25000,        0.50000,        1.00000
+316,              0.75000,        0.25000,        1.00000
+320,              0.50000,        0.25000,        1.00000
+324,              0.25000,        0.25000,        1.00000
+337,              0.75000,        0.75000,        0.50000
+341,              0.50000,        0.75000,        0.50000
+345,              0.25000,        0.75000,        0.50000
+351,              0.75000,        0.50000,        0.50000
+355,              0.50000,        0.50000,        0.50000
+359,              0.25000,        0.50000,        0.50000
+365,              0.75000,        0.25000,        0.50000
+369,              0.50000,        0.25000,        0.50000
+373,              0.25000,        0.25000,        0.50000
+438,              1.00000,        0.00000,        0.25000
+439,              1.00000,        0.00000,        0.50000
+440,              0.66667,        0.00000,        0.50000
+444,              0.66667,        0.00000,        0.25000
+447,              1.00000,        0.00000,        0.00000
+450,              0.66667,        0.00000,        0.00000
+453,              0.33333,        0.00000,        0.50000
+456,              0.33333,        0.00000,        0.25000
+460,              0.33333,        0.00000,        0.00000
+463,              0.00000,        0.00000,        0.50000
+466,              0.00000,        0.00000,        0.25000
+470,              0.00000,        0.00000,        0.00000
+473,              0.66667,        1.00000,        0.50000
+474,              1.00000,        1.00000,        0.50000
+475,              1.00000,        1.00000,        0.25000
+479,              0.66667,        1.00000,        0.25000
+482,              1.00000,        1.00000,        0.00000
+485,              0.66667,        1.00000,        0.00000
+488,              0.33333,        1.00000,        0.50000
+491,              0.33333,        1.00000,        0.25000
+495,              0.33333,        1.00000,        0.00000
+498,              0.00000,        1.00000,        0.50000
+501,              0.00000,        1.00000,        0.25000
+505,              0.00000,        1.00000,        0.00000
+508,              1.00000,        0.66667,        0.50000
+511,              1.00000,        0.66667,        0.25000
+515,              1.00000,        0.66667,        0.00000
+518,              1.00000,        0.33333,        0.50000
+521,              1.00000,        0.33333,        0.25000
+525,              1.00000,        0.33333,        0.00000
+533,              0.00000,        0.66667,        0.50000
+536,              0.00000,        0.66667,        0.25000
+540,              0.00000,        0.66667,        0.00000
+543,              0.00000,        0.33333,        0.50000
+546,              0.00000,        0.33333,        0.25000
+550,              0.00000,        0.33333,        0.00000
+559,              0.66667,        0.66667,        0.50000
+563,              0.33333,        0.66667,        0.50000
+569,              0.66667,        0.33333,        0.50000
+573,              0.33333,        0.33333,        0.50000
+584,              0.66667,        0.66667,        0.00000
+588,              0.33333,        0.66667,        0.00000
+594,              0.66667,        0.33333,        0.00000
+598,              0.33333,        0.33333,        0.00000
+608,              0.66006,        0.47128,        0.70078
+609,              0.19762,        0.64074,        0.81889
+610,              0.65771,        0.82813,        0.74829
+611,              0.59993,        0.15953,        0.76106
+612,              0.40915,        0.16311,        0.74937
+613,              0.50000,        0.50000,        0.25000
+**
+**ELSET COUNT = 172
+**HWCOLOR COMP         54     0
+*ELEMENT, TYPE=C3D4, ELSET=UPPER
+       570,       252,       262,       255,       609
+       571,       252,       310,       262,       609
+       572,       252,       296,       310,       609
+       573,       200,       252,       255,       609
+       574,       296,       306,       310,       609
+       575,       262,       310,       265,       609
+       576,       259,       265,       269,       359
+       577,       184,       337,       220,       610
+       578,       259,       345,       265,       359
+       579,       288,       608,       292,       610
+       580,       187,       288,       292,       610
+       581,       172,       288,       187,       610
+       582,       178,       172,       187,       610
+       583,       178,       220,       172,       610
+       584,       172,       220,       288,       610
+       585,       230,       240,       237,       316
+       586,       178,       184,       220,       610
+       587,       178,       190,       184,       610
+       588,       178,       187,       190,       610
+       589,       187,       292,       190,       610
+       590,       190,       292,       609,       610
+       591,       190,       609,       341,       610
+       592,       190,       341,       194,       610
+       593,       184,       190,       194,       610
+       594,       184,       194,       337,       610
+       595,       194,       341,       337,       610
+       596,       337,       341,       608,       610
+       597,       230,       337,       608,       610
+       598,       230,       608,       288,       610
+       599,       220,       230,       288,       610
+       600,       220,       337,       230,       610
+       601,       341,       609,       608,       610
+       602,       292,       608,       609,       610
+       603,       306,       316,       320,       611
+       604,       306,       608,       316,       611
+       605,       365,       608,       369,       611
+       606,       240,       608,       365,       611
+       607,       240,       316,       608,       611
+       608,       240,       133,       316,       611
+       609,       240,       365,       133,       611
+       610,       365,       139,       133,       611
+       611,       365,       369,       139,       611
+       612,       369,       149,       139,       611
+       613,       149,       145,       611,       612
+       614,       240,       127,       128,       316
+       615,       320,       129,       142,       611
+       616,       316,       129,       320,       611
+       617,       316,       133,       129,       611
+       618,       129,       133,       142,       611
+       619,       142,       133,       145,       611
+       620,       139,       145,       133,       611
+       621,       139,       149,       145,       611
+       622,       310,       320,       324,       612
+       623,       306,       320,       310,       612
+       624,       306,       611,       320,       612
+       625,       306,       608,       611,       612
+       626,       306,       609,       608,       612
+       627,       306,       310,       609,       612
+       628,       265,       609,       310,       612
+       629,       265,       310,       324,       612
+       630,       265,       324,       275,       612
+       631,       265,       275,       373,       612
+       632,       265,       373,       359,       612
+       633,       265,       359,       609,       612
+       634,       355,       609,       359,       612
+       635,       355,       608,       609,       612
+       636,       355,       369,       608,       612
+       637,       369,       611,       608,       612
+       638,       214,       255,       259,       345
+       639,       142,       611,       145,       612
+       640,       320,       611,       142,       612
+       641,       320,       142,       324,       612
+       642,       324,       142,       152,       612
+       643,       324,       152,       155,       612
+       644,       275,       324,       155,       612
+       645,       275,       155,       373,       612
+       646,       373,       155,       159,       612
+       647,       373,       159,       149,       612
+       648,       369,       373,       149,       612
+       649,       359,       373,       369,       612
+       650,       355,       359,       369,       612
+       651,       369,       149,       611,       612
+       652,       149,       155,       145,       612
+       653,       149,       159,       155,       612
+       654,       152,       145,       155,       612
+       655,       142,       145,       152,       612
+       656,       237,       240,       128,       316
+       657,       244,       136,       127,       365
+       658,       240,       365,       127,       133
+       659,       240,       351,       244,       365
+       660,       220,       224,       230,       337
+       661,       227,       288,       230,       302
+       662,       200,       207,       252,       296
+       663,       240,       127,       316,       133
+       664,       230,       337,       234,       351
+       665,       275,       162,       165,       155
+       666,       172,       220,       217,       288
+       667,       200,       252,       210,       255
+       668,       190,       292,       197,       296
+       669,       230,       234,       240,       351
+       670,       204,       210,       214,       345
+       671,       210,       255,       214,       345
+       672,       275,       165,       373,       155
+       673,       172,       174,       217,       220
+       674,       200,       204,       341,       345
+       675,       240,       244,       127,       365
+       676,       220,       230,       227,       288
+       677,       272,       275,       324,       155
+       678,       272,       152,       162,       155
+       679,       200,       255,       210,       345
+       680,       190,       200,       194,       341
+       681,       197,       207,       200,       296
+       682,       190,       197,       200,       296
+       683,       279,       165,       169,       159
+       684,       200,       207,       210,       252
+       685,       224,       234,       230,       337
+       686,       187,       197,       190,       292
+       687,       275,       165,       279,       373
+       688,       265,       275,       269,       373
+       689,       272,       324,       152,       155
+       690,       262,       265,       310,       324
+       691,       265,       272,       275,       324
+       692,       184,       224,       220,       337
+       693,       272,       162,       275,       155
+       694,       262,       272,       265,       324
+       695,       178,       174,       172,       220
+       696,       136,       365,       139,       133
+       697,       227,       230,       237,       302
+       698,       178,       184,       181,       220
+       699,       172,       174,       173,       217
+       700,       194,       200,       204,       341
+       701,       234,       244,       240,       351
+       702,       200,       210,       204,       345
+       703,       230,       237,       302,       316
+       704,       128,       316,       127,       129
+       705,       217,       220,       227,       288
+       706,       127,       365,       136,       133
+       707,       265,       269,       359,       373
+       708,       279,       373,       165,       159
+       709,       165,       159,       373,       155
+       710,       269,       275,       279,       373
+       711,       127,       129,       316,       133
+       712,       181,       220,       184,       224
+       713,       178,       181,       174,       220
+       714,       292,       302,       306,       608
+       715,       288,       302,       292,       608
+       716,       230,       302,       288,       608
+       717,       230,       316,       302,       608
+       718,       230,       240,       316,       608
+       719,       230,       351,       240,       608
+       720,       230,       337,       351,       608
+       721,       302,       316,       306,       608
+       722,       240,       351,       365,       608
+       723,       351,       355,       365,       608
+       724,       341,       355,       351,       608
+       725,       337,       341,       351,       608
+       726,       355,       369,       365,       608
+       727,       341,       355,       608,       609
+       728,       292,       608,       306,       609
+       729,       292,       306,       296,       609
+       730,       190,       292,       296,       609
+       731,       190,       200,       341,       609
+       732,       190,       296,       200,       609
+       733,       200,       296,       252,       609
+       734,       200,       255,       345,       609
+       735,       200,       345,       341,       609
+       736,       341,       345,       355,       609
+       737,       345,       359,       355,       609
+       738,       265,       359,       345,       609
+       739,       255,       265,       345,       609
+       740,       255,       265,       259,       345
+       741,       255,       262,       265,       609
+**
+**ELSET COUNT = 92
+**HWCOLOR COMP          1     0
+*ELEMENT, TYPE=C3D4, ELSET=LOWER
+       742,       456,       460,       466,       598
+       743,       444,       450,       456,       598
+       744,       444,       525,       447,       594
+       745,       536,       540,       588,       598
+       746,       453,       456,       463,       543
+       747,       533,       543,       536,       573
+       748,       444,       456,       453,       573
+       749,       466,       550,       546,       598
+       750,       536,       543,       546,       573
+       751,       491,       536,       495,       588
+       752,       438,       447,       444,       525
+       753,       444,       594,       450,       598
+       754,       444,       521,       525,       594
+       755,       438,       444,       440,       521
+       756,       540,       546,       550,       598
+       757,       453,       456,       543,       573
+       758,       488,       491,       559,       563
+       759,       505,       536,       540,       588
+       760,       440,       444,       453,       573
+       761,       491,       495,       584,       588
+       762,       444,       447,       450,       594
+       763,       485,       491,       495,       584
+       764,       440,       569,       444,       573
+       765,       495,       501,       505,       536
+       766,       473,       479,       511,       559
+       767,       438,       440,       439,       521
+       768,       491,       533,       536,       563
+       769,       444,       569,       521,       594
+       770,       473,       479,       475,       511
+       771,       491,       498,       501,       533
+       772,       495,       536,       505,       588
+       773,       438,       444,       521,       525
+       774,       473,       511,       508,       559
+       775,       440,       521,       444,       569
+       776,       511,       515,       521,       584
+       777,       491,       498,       533,       563
+       778,       515,       525,       521,       584
+       779,       460,       550,       466,       598
+       780,       439,       521,       440,       569
+       781,       456,       463,       543,       546
+       782,       479,       482,       511,       515
+       783,       533,       536,       563,       573
+       784,       479,       515,       511,       584
+       785,       450,       460,       456,       598
+       786,       521,       569,       559,       594
+       787,       559,       594,       569,       613
+       788,       559,       584,       594,       613
+       789,       521,       559,       584,       594
+       790,       536,       546,       540,       598
+       791,       511,       521,       559,       584
+       792,       479,       559,       491,       584
+       793,       491,       559,       563,       584
+       794,       511,       521,       518,       559
+       795,       491,       584,       563,       588
+       796,       563,       588,       584,       613
+       797,       559,       563,       584,       613
+       798,       559,       569,       563,       613
+       799,       563,       569,       573,       613
+       800,       536,       563,       573,       588
+       801,       563,       573,       588,       613
+       802,       573,       598,       588,       613
+       803,       491,       501,       495,       536
+       804,       536,       573,       546,       598
+       805,       536,       588,       573,       598
+       806,       588,       598,       594,       613
+       807,       584,       588,       594,       613
+       808,       444,       569,       594,       598
+       809,       521,       584,       525,       594
+       810,       444,       573,       569,       598
+       811,       569,       598,       573,       613
+       812,       569,       594,       598,       613
+       813,       444,       456,       573,       598
+       814,       439,       518,       521,       569
+       815,       456,       546,       543,       573
+       816,       479,       485,       482,       515
+       817,       508,       511,       518,       559
+       818,       456,       546,       573,       598
+       819,       488,       498,       491,       563
+       820,       479,       485,       515,       584
+       821,       491,       563,       536,       588
+       822,       473,       488,       479,       559
+       823,       479,       491,       485,       584
+       824,       491,       533,       501,       536
+       825,       479,       511,       559,       584
+       826,       479,       488,       491,       559
+       827,       475,       479,       482,       511
+       828,       518,       559,       521,       569
+       829,       473,       475,       508,       511
+       830,       473,       475,       474,       508
+       831,       460,       470,       466,       550
+       832,       456,       466,       463,       546
+       833,       456,       466,       546,       598
+**
+**ELSET COUNT = 18
+*SURFACE,TYPE=ELEMENT,NAME=LOWER_TO_UPPER
+799,S1
+798,S1
+764,S2
+814,S2
+747,S2
+828,S2
+819,S2
+822,S2
+777,S3
+746,S4
+757,S4
+760,S4
+780,S4
+783,S4
+817,S4
+758,S4
+774,S4
+830,S4
+**
+**ELSET COUNT = 18
+*SURFACE,TYPE=ELEMENT,NAME=LOWER_BOTTOM
+806,S1
+807,S1
+831,S2
+779,S2
+785,S2
+778,S2
+753,S3
+762,S3
+744,S3
+745,S3
+809,S3
+761,S3
+820,S3
+816,S3
+756,S4
+759,S4
+772,S4
+763,S4
+**
+**ELSET COUNT = 12
+*SURFACE,TYPE=ELEMENT,NAME=LOWER_SYM13
+831,S1
+742,S1
+832,S1
+746,S1
+785,S1
+743,S1
+748,S1
+760,S1
+762,S1
+752,S1
+755,S1
+767,S1
+**
+**ELSET COUNT = 12
+*SURFACE,TYPE=ELEMENT,NAME=LOWER_SYM23
+749,S1
+756,S1
+790,S1
+750,S1
+747,S1
+759,S1
+831,S3
+832,S3
+781,S3
+765,S3
+824,S3
+771,S3
+**
+**ELSET COUNT = 32
+*SURFACE,TYPE=ELEMENT,NAME=UPPER_TOP
+678,S1
+689,S1
+642,S1
+641,S1
+615,S1
+616,S1
+622,S1
+623,S1
+603,S1
+721,S1
+571,S1
+572,S1
+574,S1
+729,S1
+714,S1
+715,S1
+580,S1
+581,S1
+704,S2
+694,S2
+661,S2
+681,S2
+686,S2
+703,S3
+662,S3
+668,S3
+656,S4
+690,S4
+697,S4
+705,S4
+666,S4
+699,S4
+**
+**ELSET COUNT = 32
+*SURFACE,TYPE=ELEMENT,NAME=UPPER_TO_LOWER
+647,S1
+648,S1
+612,S1
+611,S1
+696,S1
+649,S1
+650,S1
+726,S1
+723,S1
+737,S1
+736,S1
+724,S1
+725,S1
+595,S1
+594,S1
+708,S2
+657,S2
+701,S2
+578,S2
+685,S2
+692,S2
+707,S3
+659,S3
+664,S3
+674,S3
+683,S4
+710,S4
+576,S4
+638,S4
+670,S4
+700,S4
+712,S4
+**
+**ELSET COUNT = 16
+*SURFACE,TYPE=ELEMENT,NAME=UPPER_SYM13
+653,S1
+652,S1
+654,S1
+655,S1
+621,S1
+620,S1
+619,S1
+618,S1
+709,S2
+711,S2
+683,S3
+665,S3
+678,S3
+696,S4
+706,S4
+704,S4
+**
+**ELSET COUNT = 16
+*SURFACE,TYPE=ELEMENT,NAME=UPPER_SYM23
+683,S1
+687,S1
+665,S1
+693,S1
+710,S1
+688,S1
+691,S1
+694,S1
+576,S1
+740,S1
+741,S1
+570,S1
+638,S1
+671,S1
+667,S3
+684,S3
+**
+**Property Definitions
+**
+*SOLID SECTION, ELSET=UPPER, MATERIAL=Def_Material
+*SOLID SECTION, ELSET=LOWER, MATERIAL=Def_Material
+**
+**Material Definitions
+**
+**Material:Def_Material
+*MATERIAL,NAME=Def_Material
+*ELASTIC,TYPE=ISO
+2.08000e+005,3.00000e-001
+*DENSITY
+7.80000e-009,
+*SPECIFIC HEAT
+5.00000e-001
+*CONDUCTIVITY
+4.98100e-002
+**
\ No newline at end of file