diff --git a/.gitignore b/.gitignore
index d09cfac..db3a5e9 100644
--- a/.gitignore
+++ b/.gitignore
@@ -3,6 +3,8 @@
 .ipynb_checkpoints
 docs/build/
 docs/site/
+docs/src/examples
+docs/src/packages
 *.swp
 *.lnk
 *.mess
diff --git a/docs/make.jl b/docs/make.jl
index 6769363..54288e0 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -3,24 +3,59 @@
 
 using Documenter, JuliaFEM
 
+if !haskey(Pkg.installed(), "Literate")
+    Pkg.add("Literate")
+end
+using Literate
+
+if haskey(ENV, "TRAVIS") && get(ENV, "INSTALL_PYPLOT", "false") == "true"
+    info("inside TRAVIS, installing PyPlot + matplotlib")
+    Pkg.add("PyPlot")
+    run(`pip install matplotlib`)
+end
+
 """
     copy_docs(pkg_name)
 
 Copy documentation of some package `pkg_name` to `docs/src/pkg_name`,
 where `pkg_name` is the name of package. Return true if success, false
 otherwise.
+
+If package is undocumented, i.e. directory `docs/src` is missing,
+but there still exists `README.md`, copy that file to
+`docs/src/pkg_name/index.md`.
 """
 function copy_docs(pkg_name)
-    src = joinpath(Pkg.dir(pkg_name), "docs", "src")
-    if !isdir(src)
-        warn("Cannot copy documentation of package $pkg_name from $src: ",
-             "No such directory exists. (Is the package in REQUIRE of JuliaFEM?)")
-        return false
+
+    src_dir = Pkg.dir(pkg_name, "docs", "src")
+    pkg_dir = Pkg.dir("JuliaFEM", "docs", "src", "packages")
+    dst_dir = joinpath(pkg_dir, pkg_name)
+    isdir(pkg_dir) || mkpath(pkg_dir)
+
+    # if can find pkg_name/docs/src =>
+    # copy that to docs/src/packages/pkg_name
+    if isdir(src_dir)
+        isdir(dst_dir) || cp(src_dir, dst_dir)
+        info("Copied documentation of package $pkg_name from $src_dir succesfully.")
+        return true
     end
-    dst = joinpath(Pkg.dir("JuliaFEM"), "docs", "src", pkg_name)
-    cp(src, dst; remove_destination=true)
-    info("Copied documentation of package $pkg_name from $src succesfully.")
-    return true
+
+    # if can find pkg_name/README.md =>
+    # copy that to docs/src/packages/pkg_name/index.md
+    readme_file = Pkg.dir(pkg_name, "README.md")
+    if isfile(readme_file)
+        isdir(dst_dir) || mkpath(dst_dir)
+        dst_file = joinpath(dst_dir, "index.md")
+        isfile(dst_file) || cp(readme_file, dst_file)
+        info("Copied README.md of package $pkg_name from $readme_file succesfully.")
+        return true
+    end
+
+    warn("Cannot copy documentation of package $pkg_name from $src_dir: ",
+         "No such directory exists. (Is the package in REQUIRE of JuliaFEM?)")
+
+    return false
+
 end
 
 """
@@ -48,22 +83,19 @@ add_page!(PACKAGES, "Theory" => "MyPackage/theory.md")
 """
 function add_page!(dst, src)
     file = isa(src, Pair) ? src.second : src
-    if !isfile(joinpath(Pkg.dir("JuliaFEM"), "docs", "src", file))
-        warn("Cannot add page $file: no such file")
-        return false
+    src_dir = Pkg.dir("JuliaFEM", "docs", "src")
+    if isfile(joinpath(src_dir, file))
+        push!(dst, src)
+        return true
     end
-    push!(dst, src)
-    return true
+    if isfile(joinpath(src_dir, "packages", file))
+        push!(dst, joinpath("packages", src))
+        return true
+    end
+    warn("Cannot add page $file: no such file")
+    return false
 end
 
-#=
-if haskey(ENV, "TRAVIS")
-    println("inside TRAVIS, installing PyPlot + matplotlib")
-    Pkg.add("PyPlot")
-    run(`pip install matplotlib`)
-end
-=#
-
 USER_GUIDE = []
 
 # Developer's guide is published in FEMBase.jl
@@ -82,18 +114,61 @@ if copy_docs("FEMBase")
 end
 
 # Let's construct here some description for packages
+
 PACKAGES = []
+copy_docs("FEMBase")                && add_page!(PACKAGES, "FEMBase/index.md")
+copy_docs("FEMBasis")               && add_page!(PACKAGES, "FEMBasis/index.md")
+copy_docs("FEMQuad")                && add_page!(PACKAGES, "FEMQuad/index.md")
+copy_docs("FEMSparse")              && add_page!(PACKAGES, "FEMSparse/index.md")
+copy_docs("Materials")              && add_page!(PACKAGES, "Materials/index.md")
+copy_docs("AsterReader")            && add_page!(PACKAGES, "AsterReader/index.md")
+copy_docs("AbaqusReader")           && add_page!(PACKAGES, "AbaqusReader/index.md")
+copy_docs("LinearImplicitDynamics") && add_page!(PACKAGES, "LinearImplicitDynamics/index.md")
+copy_docs("HeatTransfer")           && add_page!(PACKAGES, "HeatTransfer/index.md")
+copy_docs("PlaneElasticity")        && add_page!(PACKAGES, "PlaneElasticity/index.md")
+copy_docs("FEMBeam")                && add_page!(PACKAGES, "FEMBeam/index.md")
+copy_docs("FEMCoupling")            && add_page!(PACKAGES, "FEMCoupling/index.md")
+copy_docs("FEMTruss")               && add_page!(PACKAGES, "FEMTruss/index.md")
+copy_docs("Mortar2D")               && add_page!(PACKAGES, "Mortar2D/index.md")
+copy_docs("Mortar3D")               && add_page!(PACKAGES, "Mortar3D/index.md")
+copy_docs("MortarContact2D")        && add_page!(PACKAGES, "MortarContact2D/index.md")
+copy_docs("MortarContact2DAD")      && add_page!(PACKAGES, "MortarContact2DAD/index.md")
+copy_docs("OptoMechanics")          && add_page!(PACKAGES, "OptoMechanics/index.md")
+copy_docs("Miniball")               && add_page!(PACKAGES, "Miniball/index.md")
+copy_docs("ModelReduction")         && add_page!(PACKAGES, "ModelReduction/index.md")
+copy_docs("NodeNumbering")          && add_page!(PACKAGES, "NodeNumbering/index.md")
+#copy_docs("Xdmf")                   && add_page!(PACKAGES, "Xdmf/index.md")
+copy_docs("UMAT")                   && add_page!(PACKAGES, "UMAT/index.md")
 
-if copy_docs("FEMQuad")
-    add_page!(PACKAGES, "FEMQuad/index.md")
+# Generate examples using Literate.jl
+
+function license_stripper(s)
+    # return replace(s, r"# This file is a part of JuliaFEM.*?\n\n"sm, "")
+    lines = split(s, '\n')
+    function myfilt(line)
+        contains(line, "# This file is a part of JuliaFEM.") && return false
+        contains(line, "# License is MIT") && return false
+        return true
+    end
+    new_s = join(filter(myfilt, lines), '\n')
+    return new_s
 end
-
-if copy_docs("FEMBasis")
-    add_page!(PACKAGES, "FEMBasis/index.md")
-end
-
-if copy_docs("HeatTransfer")
-    add_page!(PACKAGES, "HeatTransfer/index.md")
+EXAMPLES = []
+docs_src = Pkg.dir("JuliaFEM", "docs", "src")
+ex_src = Pkg.dir("JuliaFEM", "examples")
+ex_dst = joinpath(docs_src, "examples")
+for ex_file in readdir(ex_src)
+    endswith(ex_file, ".jl") || continue
+    startswith(ex_file, "test_") && continue
+    dst_file = Literate.markdown(joinpath(ex_src, ex_file), ex_dst;
+                                 documenter=true, preprocess=license_stripper)
+    ex_dir = joinpath(ex_src, first(splitext(ex_file)))
+    dst_ex_dir = joinpath(ex_dst, first(splitext(ex_file)))
+    if isdir(ex_dir) && !isdir(dst_ex_dir)
+        info("Copy $ex_dir to $dst_ex_dir")
+        cp(ex_dir, dst_ex_dir)
+    end
+    add_page!(EXAMPLES, relpath(dst_file, docs_src))
 end
 
 # API documentation
@@ -103,6 +178,7 @@ LIBRARY = ["api.md"]
 PAGES = []
 push!(PAGES, "Home" => "index.md")
 #push!(PAGES, "User's guide" => USER_GUIDE)
+push!(PAGES, "Examples" => EXAMPLES)
 push!(PAGES, "Developer's guide" => DEVELOPER_GUIDE)
 push!(PAGES, "Description of packages" => PACKAGES)
 #push!("API documentation" => LIBRARY)
diff --git a/examples/2d_hertz_contact.jl b/examples/2d_hertz_contact.jl
new file mode 100644
index 0000000..4dcbe7f
--- /dev/null
+++ b/examples/2d_hertz_contact.jl
@@ -0,0 +1,221 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
+
+# # 2D Hertz contact problem
+
+# ![](2d_hertz_contact/model.png)
+
+# In the example, a cylinder is pressed agains block with a force of 35 kN.
+# A similar example can be found from NAFEMS report FENET D3613 (advanced
+# finite element contact benchmarks).
+#
+# Solution for maximum pressure ``p_0`` and contact radius ``a`` is
+# ```math
+#   p_{0} = \sqrt{\frac{FE}{2\pi R}}, \\
+#   a     = \sqrt{\frac{8FR}{\pi E}},
+# ```
+# where
+#
+# ```math
+#   E = \frac{2E_{1}E_{2}}{E_{2}\left(1-\nu_{1}^{2}\right)+E_{1}\left(1-\nu_{2}^{2}\right)}.
+# ```
+#
+# Substituting values, one gets accurate solution to be ``p_0 = 3585 \;\mathrm{MPa}`` and
+# ``a = 6.21 \;\mathrm{mm}``.
+
+using JuliaFEM
+using JuliaFEM.Preprocess
+using JuliaFEM.Postprocess
+using Logging
+Logging.configure(level=INFO)
+add_elements! = JuliaFEM.add_elements!
+
+# Simulation starts by reading the mesh. Model is constructed and meshed using
+# SALOME, thus mesh format is .med. Mesh type is quite simple structure,
+# containing things like `mesh.nodes`, `mesh.elements` and so on. Keep on mind,
+# that Mesh contains only standard Julia types and we think it as a structure
+# helping us to construct elements needed in simulation. In principle, we don't
+# need to use `Mesh` in simulation anyway if we figure some other way to define
+# the geometry for elements.
+
+datadir = Pkg.dir("JuliaFEM", "examples", "2d_hertz_contact")
+meshfile = joinpath(datadir, "hertz_2d_full.med")
+mesh = aster_read_mesh(meshfile)
+for (elset_name, element_ids) in mesh.element_sets
+    nel = length(element_ids)
+    println("Element set $elset_name contains $nel elements.")
+end
+for (nset_name, node_ids) in mesh.node_sets
+    nno = length(node_ids)
+    println("Node set $nset_name contains $nno nodes.")
+end
+nnodes = length(mesh.nodes)
+println("Total number of nodes in mesh: $nnodes")
+nelements = length(mesh.elements)
+println("Total number of elements in mesh: $nelements")
+
+# Next, define two bodies. Technically, we could have only one problem and add
+# elements from both bodies to the same problem, but defining two different
+# problems is recommended for clarity. Plain strain assumption is used.
+# To make clear what is happening here: we first create a set of elements
+# (elements are in vector called `upper_elements`), then we define new
+# problem which type is `Elasticity`, give it some meaningful name (this time
+# `cylinder`), and last value 2 means that problems does have two degrees of
+# freedom per node.
+
+upper_elements = create_elements(mesh, "CYLINDER")
+update!(upper_elements, "youngs modulus", 70.0e3)
+update!(upper_elements, "poissons ratio", 0.3)
+upper = Problem(Elasticity, "cylinder", 2)
+upper.properties.formulation = :plane_strain
+add_elements!(upper, upper_elements)
+
+lower_elements = create_elements(mesh, "BLOCK")
+update!(lower_elements, "youngs modulus", 210.0e3)
+update!(lower_elements, "poissons ratio", 0.3)
+lower = Problem(Elasticity, "block", 2)
+lower.properties.formulation = :plane_strain
+add_elements!(lower, lower_elements)
+
+# Next we define some boundary conditions: creating "boundary" problems goes
+# in the same way than defining "field" problems, the only difference is that
+# we add extra argument giving what field are we tring to fix. This time,
+# we have 2 dofs / node and we fix displacement in direction 2.
+
+bc_fixed_elements = create_elements(mesh, "FIXED")
+update!(bc_fixed_elements, "displacement 2", 0.0)
+bc_fixed = Problem(Dirichlet, "fixed", 2, "displacement")
+add_elements!(bc_fixed, bc_fixed_elements)
+
+# Defining symmetry boundary condition goes with the same idea
+
+bc_sym_23_elements = create_elements(mesh, "SYM23")
+update!(bc_sym_23_elements, "displacement 1", 0.0)
+bc_sym_23 = Problem(Dirichlet, "symmetry line 23", 2, "displacement")
+add_elements!(bc_sym_23, bc_sym_23_elements)
+
+# Next we define point load. To define that, we first need to find some node
+# near the top of cylinder, using function `find_nearest_node`. Then we create
+# a new problem, again of type Elasticity. Like told already, we don't need to
+# use `Mesh` if we have some other procedure to define the geometry of the
+# element (and it's connectivity, of course). So we can directly create an
+# element of type `Poi1`, meaning 1-node point element, update it's geometry
+# and apply 35.0e3 kN load in negative y-direction:
+
+nid = find_nearest_node(mesh, [0.0, 100.0])
+load = Problem(Elasticity, "point load", 2)
+load.properties.formulation = :plane_strain
+load.elements = [Element(Poi1, [nid])]
+update!(load.elements, "geometry", mesh.nodes)
+update!(load.elements, "displacement traction force 2", -35.0e3)
+
+# Next, we define another boudary problem, this time the type of problem is
+# Contact2D, which is a mortar contact formulation for two dimensions.
+# Elements are added using `add_slave_elements!` and `add_master_elements!`.
+# Problems, in general, can have some properties defined, like the formulation
+# in `Elasticity` (we also have `:plane_stress`). For contact, we need to swap
+# normal direction for meshes created by SALOME because in Code Aster, element
+# orientation is defined opposite to what is used in ABAQUS, and in JuliaFEM in
+# general we follow the same conventions what are used in ABAQUS.
+
+contact = Problem(Contact2D, "contact", 2, "displacement")
+contact.properties.rotate_normals = true
+contact_slave_elements = create_elements(mesh, "BLOCK_TO_CYLINDER")
+contact_master_elements = create_elements(mesh, "CYLINDER_TO_BLOCK")
+add_master_elements!(contact, contact_master_elements)
+add_slave_elements!(contact, contact_slave_elements)
+
+# After all problems are defined, we define some `Analysis`, which can be e.g.
+# static analysis, dynamic analysis, modal analysis, linear perturbation
+# analysis and so on. Here, the analysis type is `Nonlinear`, which is nonlinear
+# quasistatic analysis. In the same manner as we do `add_elements!` to add
+# elements to `Problem`, we use `add_problems!` to add problems to analysis.
+# Because we are not restricted to some particular input and output formats,
+# we "connect" a `ResultsWriter` to our analysis, this time we want to visualize
+# results using ParaView, thus we write our results to Xdmf format, which uses
+# well defined standards XML and HDF to store model data.
+
+step = Analysis(Nonlinear)
+add_problems!(step, [upper, lower, bc_fixed, bc_sym_23, load, contact])
+xdmf = Xdmf("2d_hertz_results"; overwrite=true)
+## todo for 2d
+## for body in (upper, lower)
+##     push!(body.postprocess_fields, "stress")
+## end
+add_results_writer!(step, xdmf)
+
+# In last part, we run the analysis.
+
+step()
+
+# # Results
+
+# Results are stored in `2d_hertz_results.xmf` and `2d_hertz_results.h5` for
+# visual inspection. We can also postprocess results programmatically because
+# we are inside a real scripting / programming environment all the time.  For
+# example, we can integrate the resultant force in normal and tangential direction
+# in contact surface to validate our result.
+
+Rn = 0.0
+Rt = 0.0
+time = 0.0
+for sel in contact_slave_elements
+    for ip in get_integration_points(sel)
+        w = ip.weight*sel(ip, time, Val{:detJ})
+        n = sel("normal", ip, time)
+        t = sel("tangent", ip, time)
+        la = sel("lambda", ip, time)
+        Rn += w*dot(n, la)
+        Rt += w*dot(t, la)
+    end
+end
+
+println("2d hertz contact resultant forces: Rn = $Rn, Rt = $Rt")
+
+using Base.Test
+@test isapprox(Rn, 35.0e3)
+@test isapprox(Rt, 0.0)
+
+# Visualization of the results can be done using ParaView:
+# ![](2d_hertz_contact/results_displacement.png)
+
+# For optimization loops, we want to programmatically find, for example, maximum
+# contact pressure. We can, for example, get all the values in nodes:
+
+lambda = contact("lambda", time)
+normal = contact("normal", time)
+p0 = 0.0
+p0_acc = 3585.0
+for (nid, n) in normal
+    lan = dot(n, lambda[nid])
+    println("$nid => $lan")
+    p0 = max(p0, lan)
+end
+p0 = round(p0, 2)
+rtol = round(norm(p0-p0_acc)/max(p0,p0_acc)*100, 2)
+println("Maximum contact pressure p0 = $p0, p0_acc = $p0_acc, rtol = $rtol %")
+
+# To get rough approximation where does the contact open, we can find the element
+# from slave contact surface, where contact pressure is zero in the other node
+# and something nonzero in the other node.
+
+a_rad = 0.0
+for element in contact_slave_elements
+    la1, la2 = element("lambda", time)
+    p1, p2 = norm(la1), norm(la2)
+    a, b = isapprox(p1, 0.0), isapprox(p2, 0.0)
+    if (a && !b) || (b && !a)
+        X1, X2 = element("geometry", time)
+        println("Contact opening element geometry: X1 = $X1, X2 = $X2")
+        println("Contact opening element lambda: la1 = $la1, la2 = $la2")
+        x11, y11 = X1
+        x12, y12 = X2
+        a_rad = 1/2*abs(x11+x12)
+        break
+    end
+end
+println("Contact radius: $a_rad")
+
+# This example briefly described some of the core features of JuliaFEM.
+
+close(xdmf.hdf) # src
diff --git a/test/test_problems_contact_2d/hertz_2d_full.med b/examples/2d_hertz_contact/hertz_2d_full.med
similarity index 100%
rename from test/test_problems_contact_2d/hertz_2d_full.med
rename to examples/2d_hertz_contact/hertz_2d_full.med
diff --git a/examples/2d_hertz_contact/model.png b/examples/2d_hertz_contact/model.png
new file mode 100644
index 0000000..e7ecd68
Binary files /dev/null and b/examples/2d_hertz_contact/model.png differ
diff --git a/examples/2d_hertz_contact/results_displacement.png b/examples/2d_hertz_contact/results_displacement.png
new file mode 100644
index 0000000..6c6f6f7
Binary files /dev/null and b/examples/2d_hertz_contact/results_displacement.png differ
diff --git a/examples/generate_stiffness_matrices.jl b/examples/generate_stiffness_matrices.jl
new file mode 100644
index 0000000..65497e9
--- /dev/null
+++ b/examples/generate_stiffness_matrices.jl
@@ -0,0 +1,19 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
+
+# # Generating local matrices for problems
+
+using JuliaFEM
+
+# Plane stress Quad4 element with linear material model:
+
+X = Dict(1 => [0.0, 0.0], 2 => [2.0, 0.0], 3 => [2.0, 2.0], 4 => [0.0, 2.0])
+element = Element(Quad4, [1, 2, 3, 4])
+update!(element, "geometry", X)
+update!(element, "youngs modulus", 288.0)
+update!(element, "poissons ratio", 1/3)
+problem = Problem(Elasticity, "test problem", 2)
+problem.properties.formulation = :plane_stress
+add_elements!(problem, [element])
+assemble!(problem, 0.0)
+K = round.(full(problem.assembly.K), 5)
diff --git a/test/test_problems_contact_2d.jl b/test/test_problems_contact_2d.jl
deleted file mode 100644
index f17f14c..0000000
--- a/test/test_problems_contact_2d.jl
+++ /dev/null
@@ -1,175 +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
-
-testdir = joinpath(Pkg.dir("JuliaFEM"), "test")
-datadir = first(splitext(basename(@__FILE__)))
-
-@testset "hertz contact, full 2d model, linear elements, flat slave surface" begin
-    meshfile = joinpath(testdir, 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(Contact2D, "contact between block and cylinder", 2, "displacement")
-    contact.properties.rotate_normals = true
-    contact_slave_elements = create_elements(mesh, "BLOCK_TO_CYLINDER")
-    contact_master_elements = create_elements(mesh, "CYLINDER_TO_BLOCK")
-    add_master_elements!(contact, contact_master_elements)
-    add_slave_elements!(contact, contact_slave_elements)
-
-    solver = Solver(Nonlinear)
-    push!(solver, upper, lower, bc_fixed, bc_sym_23, load, contact)
-    solver()
-
-    node_ids, la = get_nodal_vector(contact_slave_elements, "lambda", 0.0)
-    node_ids, n = get_nodal_vector(contact_slave_elements, "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 contact_slave_elements
-        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("lambda", 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 resultant force
-    @test isapprox(Rn, 35.0e3; rtol=0.020)
-    @test isapprox(Rt, 0.0; atol=200.0)
-end
-
-function get_model()
-    meshfile = joinpath(testdir, 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(Contact2D, "interface", 2, "displacement")
-    interface.properties.rotate_normals = true
-    interface_slave_elements = create_elements(mesh, "LOWER_TOP")
-    interface_master_elements = create_elements(mesh, "UPPER_BOTTOM")
-    add_master_elements!(interface, interface_master_elements)
-    add_slave_elements!(interface, interface_slave_elements)
-
-    # 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, lambda = get_nodal_vector(interface.elements, "lambda", 0.0)
-    u2 = [u[2] for u in displacement]
-    f2 = [f[2] for f in lambda]
-    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, lambda = get_nodal_vector(interface.elements, "lambda", 0.0)
-    u2 = [u[2] for u in displacement]
-    f2 = [f[2] for f in lambda]
-    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
deleted file mode 100644
index 55e8b7a..0000000
Binary files a/test/test_problems_contact_2d/block_2d.med and /dev/null differ
diff --git a/test/test_problems_contact_2d_hertz_1.jl b/test/test_problems_contact_2d_hertz_1.jl
deleted file mode 100644
index b9eb217..0000000
--- a/test/test_problems_contact_2d_hertz_1.jl
+++ /dev/null
@@ -1,86 +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
-
-pkgdir = Pkg.dir("JuliaFEM")
-datadir = joinpath(pkgdir, "test", 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
-
-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(Contact2D, "contact between block and cylinder", 2, "displacement")
-contact.properties.rotate_normals = true
-contact_slave_elements = create_elements(mesh, "CYLINDER_TO_BLOCK")
-contact_master_elements = create_elements(mesh, "BLOCK_TO_CYLINDER")
-add_master_elements!(contact, contact_master_elements)
-add_slave_elements!(contact, contact_slave_elements)
-
-solver = Solver(Nonlinear)
-push!(solver, upper, lower, bc_fixed, bc_sym_23, load, contact)
-solver()
-
-node_ids, la = get_nodal_vector(contact_slave_elements, "lambda", 0.0)
-node_ids, n = get_nodal_vector(contact_slave_elements, "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 contact_slave_elements
-    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("lambda", 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 resultant force
-@test isapprox(Rn, 35.0e3; rtol=0.020)
-@test isapprox(Rt, 0.0; atol=200.0)
diff --git a/test/test_problems_contact_2d_hertz_1/hertz_2d_full.med b/test/test_problems_contact_2d_hertz_1/hertz_2d_full.med
deleted file mode 100644
index 2733f0b..0000000
Binary files a/test/test_problems_contact_2d_hertz_1/hertz_2d_full.med and /dev/null differ