diff --git a/src/JuliaFEM.jl b/src/JuliaFEM.jl
index f33bd2c..3a9920a 100644
--- a/src/JuliaFEM.jl
+++ b/src/JuliaFEM.jl
@@ -6,6 +6,8 @@ This is JuliaFEM -- Finite Element Package
 """
 module JuliaFEM
 
+using Compat
+import Compat.String
 importall Base
 
 include("fields.jl")
@@ -15,10 +17,10 @@ export AbstractPoint, Point, IntegrationPoint, IP, Node
 
 ### ELEMENTS ###
 include("elements.jl") # common element routines
-export Node, AbstractElement, Element, update!, get_connectivity, get_basis, get_dbasis
+export Node, AbstractElement, Element, update!, get_connectivity, get_basis, get_dbasis, inside, get_local_coordinates
 include("elements_lagrange_macro.jl") # Continuous Galerkin (Lagrange) elements generated using macro
 include("elements_lagrange.jl") # Continuous Galerkin (Lagrange) elements
-export get_reference_coordinates
+export get_reference_coordinates, get_interpolation_polynomial
 export Poi1,
        Seg2, Seg3,
        Tri3, Tri6, Quad4, Quad8, Quad9,
diff --git a/src/preprocess.jl b/src/preprocess.jl
index 2c19c79..f7b19ab 100644
--- a/src/preprocess.jl
+++ b/src/preprocess.jl
@@ -89,13 +89,8 @@ function filter_by_element_set(mesh::Mesh, set_name::String)
 end
 
 function create_elements(mesh::Mesh)
-    elements = Element[]
-    for (elid, elcon) in mesh.elements
-        eltype = mesh.element_types[elid]
-        element = Element(JuliaFEM.(eltype), elcon)
-        update!(element, "geometry", mesh.nodes)
-        push!(elements, element)
-    end
+    elements = [Element(JuliaFEM.(mesh.element_types[elid]), elcon) for (elid, elcon) in mesh.elements]
+    update!(elements, "geometry", mesh.nodes)
     return elements
 end
 
@@ -103,7 +98,7 @@ function create_elements(mesh::Mesh, element_sets::String...)
     elements = Element[]
     for element_set in element_sets
         new_elements = create_elements(filter_by_element_set(mesh, element_set))
-        push!(elements, new_elements...)
+        elements = [elements; new_elements]
     end
     return elements
 end
diff --git a/src/preprocess_aster_reader.jl b/src/preprocess_aster_reader.jl
index 83a6bfd..67088b8 100644
--- a/src/preprocess_aster_reader.jl
+++ b/src/preprocess_aster_reader.jl
@@ -346,11 +346,13 @@ Dict containing fields "nodes" and "connectivity".
 """
 function parse_aster_med_file(fn::String, mesh_name=nothing; debug=false)
     med = MEDFile(fn)
-    if isa(mesh_name, Void)
-        mesh_names = get_mesh_names(med::MEDFile)
-        all_meshes = join(mesh_names, ", ")
+    mesh_names = get_mesh_names(med::MEDFile)
+    all_meshes = join(mesh_names, ", ")
+    if mesh_name == nothing
         length(mesh_names) == 1 || error("several meshes found from med, pick one: $all_meshes")
         mesh_name = mesh_names[1]
+    else
+        mesh_name in mesh_names || error("Mesh $mesh_name not found from mesh file $fn. Available meshes: $all_meshes")
     end
     nsets = get_node_sets(med, mesh_name)
     elsets = get_element_sets(med, mesh_name)
diff --git a/src/problems_contact_2d_autodiff.jl b/src/problems_contact_2d_autodiff.jl
index 68b1608..621c6a2 100644
--- a/src/problems_contact_2d_autodiff.jl
+++ b/src/problems_contact_2d_autodiff.jl
@@ -21,7 +21,7 @@ xi
 """
 function project_from_master_to_slave{E<:MortarElements2D}(
     slave_element::Element{E}, x1_::DVTI, n1_::DVTI, x2::Vector;
-    tol=1.0e-10, max_iterations=20)
+    tol=1.0e-10, max_iterations=20, debug=false)
 
     x1(xi1) = vec(get_basis(slave_element, [xi1], time))*x1_
     dx1(xi1) = vec(get_dbasis(slave_element, [xi1], time))*x1_
@@ -32,21 +32,32 @@ function project_from_master_to_slave{E<:MortarElements2D}(
     dR(xi1) = cross2(dx1(xi1), n1(xi1)) + cross2(x1(xi1)-x2, dn1(xi1))
 
     xi1 = 0.0
+    xi1_next = 0.0
     dxi1 = 0.0
     for i=1:max_iterations
         dxi1 = -R(xi1)/dR(xi1)
-        xi1 += dxi1
-        if norm(dxi1) < tol
-            return xi1
+        dxi1 = clamp(dxi1, -0.3, 0.3)
+        xi1_next = clamp(xi1 + dxi1, -1.0, 1.0)
+        if norm(xi1_next - xi1) < tol
+            return xi1_next
         end
+        if debug
+            info("xi1 = $xi1")
+            info("R(xi1) = $(R(xi1))")
+            info("dR(xi1) = $(dR(xi1))")
+            info("dxi1 = $dxi1")
+            info("norm = $(norm(xi1_next - xi1))")
+            info("xi1_next = $xi1_next")
+        end
+        xi1 = xi1_next
     end
 
-    info("x1 = $(ForwardDiff.get_value(x1_.data))")
-    info("n1 = $(ForwardDiff.get_value(n1_.data))")
-    info("x2 = $(ForwardDiff.get_value(x2))")
-    info("xi1 = $(ForwardDiff.get_value(xi1)), dxi1 = $(ForwardDiff.get_value(dxi1))")
-    info("-R(xi1) = $(ForwardDiff.get_value(-R(xi1)))")
-    info("dR(xi1) = $(ForwardDiff.get_value(dR(xi1)))")
+    info("x1 = $x1_")
+    info("n1 = $n1_")
+    info("x2 = $x2")
+    info("xi1 = $xi1, dxi1 = $dxi1")
+    info("-R(xi1) = $(-R(xi1))")
+    info("dR(xi1) = $(dR(xi1))")
     error("find projection from master to slave: did not converge")
 
 end
@@ -157,17 +168,8 @@ function assemble!(problem::Problem{Contact}, time::Float64,
                 #distance > props.maximum_distance && continue
 
                 # calculate segmentation: we care only about endpoints
-                # note: these are quadratic/cubic functions, analytical solution possible
-                xi1a = -Inf
-                xi1b = -Inf
-                try
-                    xi1a = project_from_master_to_slave(slave_element, x1, n1, x2[1])
-                    xi1b = project_from_master_to_slave(slave_element, x1, n1, x2[2])
-                catch
-                    info("failed to create projection!!!!")
-                    # TODO
-                    continue
-                end
+                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
diff --git a/test/test_contact_2d_finite_sliding.jl b/test/test_contact_2d_finite_sliding.jl
index 8e35a7f..b491af3 100644
--- a/test/test_contact_2d_finite_sliding.jl
+++ b/test/test_contact_2d_finite_sliding.jl
@@ -54,3 +54,18 @@ using JuliaFEM.Test
 #   @test isapprox(normu, 0.49745873784105105)
     @test isapprox(normu, 0.49745872893844145)
 end
+
+#= TODO: Fix test, this is not converging
+@testset "project from master to slave" begin
+    el = Element(Seg2, [1, 2])
+    x1 = DVTI(Vector{Float64}[
+             [ 0.07406987526791842, 0.6628967239474994],
+             [-0.24633092752656838, 0.4732606367688589]])
+    n1 = DVTI(Vector{Float64}[
+             [0.40398625635635355, 0.9147650543583191],
+             [-0.5093430176405901, 0.860563588807229]])
+    x2 = [0.5049198709043257, 0.27765317280577695]
+    xi = project_from_master_to_slave(el, x1, n1, x2; debug=true)
+    info("xi = $xi")
+end
+=#
diff --git a/test/test_elements_2.jl b/test/test_elements_2.jl
new file mode 100644
index 0000000..3bbb90c
--- /dev/null
+++ b/test/test_elements_2.jl
@@ -0,0 +1,50 @@
+using JuliaFEM
+using JuliaFEM.Test
+
+@testset "inverse isoparametric mapping" begin
+    el = Element(Quad4, [1, 2, 3, 4])
+    X = Dict{Int64, Vector{Float64}}(
+        1 => [0.0, 0.0],
+        2 => [1.0, 0.0],
+        3 => [1.0, 1.0],
+        4 => [0.0, 1.0])
+    update!(el, "geometry", X)
+    time = 0.0
+    X1 = el("geometry", [0.1, 0.2], time)
+    xi = get_local_coordinates(el, X1, time)
+    X2 = el("geometry", xi, time)
+    info("X1 = $X1, X2 = $X2")
+    @test isapprox(X1, X2)
+end
+
+@testset "inside of linear element" begin
+    el = Element(Quad4, [1, 2, 3, 4])
+    X = Dict{Int64, Vector{Float64}}(
+        1 => [0.0, 0.0],
+        2 => [1.0, 0.0],
+        3 => [1.0, 1.0],
+        4 => [0.0, 1.0])
+    update!(el, "geometry", X)
+    time = 0.0
+    @test inside(el, [0.5, 0.5], time) == true
+    @test inside(el, [1.0, 0.5], time) == true
+    @test inside(el, [1.0, 1.0], time) == true
+    @test inside(el, [1.01, 1.0], time) == false
+    @test inside(el, [1.0, 1.01], time) == false
+end
+
+@testset "inside of quadratic element" begin
+    el = Element(Tri6, [1, 2, 3, 4, 5, 6])
+    X = Dict{Int64, Vector{Float64}}(
+        1 => [0.0, 0.0],
+        2 => [1.0, 0.0],
+        3 => [0.0, 1.0],
+        4 => [0.5, 0.2],
+        5 => [0.8, 0.6],
+        6 => [-0.2, 0.5])
+    update!(el, "geometry", X)
+    p = [0.94, 0.3] # visually checked to be inside
+    @test inside(el, p, 0.0) == true
+    p = [-0.2, 0.8] # visually checked to be outside
+    @test inside(el, p, 0.0) == false
+end
diff --git a/test/test_heat.jl b/test/test_heat.jl
index db47281..02467ff 100644
--- a/test/test_heat.jl
+++ b/test/test_heat.jl
@@ -9,7 +9,7 @@ using JuliaFEM.Postprocess
 @testset "Tet10 + convection" begin
     # For some reason Tet10 fails, maybe because of convection.
     mesh_file = Pkg.dir("JuliaFEM") * "/test/testdata/primitives.med"
-    mesh = aster_read_mesh(mesh_file, "Tet10")
+    mesh = aster_read_mesh(mesh_file, "TETRA_TET10_1")
     prob = Problem(Heat, "tet", 1)
     face = Problem(Heat, "face 4", 1)
     fixed = Problem(Dirichlet, "fixed face 3", 1, "temperature")
diff --git a/test/test_heat_2.jl b/test/test_heat_2.jl
new file mode 100644
index 0000000..fb8a9cd
--- /dev/null
+++ b/test/test_heat_2.jl
@@ -0,0 +1,49 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
+using JuliaFEM
+using JuliaFEM.Preprocess
+using JuliaFEM.Postprocess
+using JuliaFEM.Test
+
+@testset "3d rod" begin
+    mesh = aster_read_mesh(Pkg.dir("JuliaFEM")*"/test/testdata/primitives.med", "CYLINDER_20_TET4")
+    problem = Problem(Heat, "rod of length 20", 1)
+    problem.elements = create_elements(mesh, "CYLINDER")
+    update!(problem, "temperature thermal conductivity", 200.0)
+    outer = Problem(Heat, "outer surface", 1)
+    outer.elements = create_elements(mesh, "FACE2", "FACE3")
+    update!(outer, "temperature external temperature", 20.0)
+    update!(outer, "temperature heat transfer coefficient", 1.0)
+    #midline = Problem(Heat, "midline of rod", 1)
+    #midline.elements = create_elements(mesh, "INNER_LINE")
+    boundary = Problem(Dirichlet, "homogeneous dirichlet boundary", 1, "temperature")
+    boundary.elements = create_elements(mesh, "FACE1")
+    update!(boundary, "temperature 1", 100.0)
+    #solver = LinearSolver(problem, outer, boundary, midline)
+    solver = LinearSolver(problem, outer, boundary)
+    solver()
+
+    L = 20
+    k = 200.0
+    Tu = 20.0
+    h = 1.0
+    P = 2*pi
+    A = pi
+    α = h
+    β = sqrt((h*P)/(k*A))
+    T0 = 100.0
+    C = [1.0 1.0; (α+k*β)*exp(β*L) (α-k*β)*exp(-β*L)] \ [T0-Tu, 0]
+    T(x) = dot(C, [exp(β*x), exp(-β*x)]) + Tu
+
+    T_diff = []
+    for x in linspace(0, 20)
+        T_FEM = problem("temperature", [x, 0.0, 0.0])[1]
+        T_ACC = T(x)
+        push!(T_diff, norm(T_FEM - T_ACC))
+        info("x = $x, T_FEM = $T_FEM, T_ACC = $T_ACC")
+    end
+    info("mean diff = ", mean(T_diff))
+    # mean diff = 1.14
+    @test mean(T_diff) < 1.2
+end
+
diff --git a/test/test_heat_3.jl b/test/test_heat_3.jl
new file mode 100644
index 0000000..0d77cf7
--- /dev/null
+++ b/test/test_heat_3.jl
@@ -0,0 +1,52 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
+using JuliaFEM
+using JuliaFEM.Preprocess
+using JuliaFEM.Postprocess
+using JuliaFEM.Test
+
+@testset "2d poisson problem with known analytical solution" begin
+    # from FENiCS tutorial, u(x,y) = 1 + x² + 2y² on [0x1]×[0,1]
+    # and u₀(x,y) = 1 + x² + 2y², f(x,y) = -6
+
+    mesh_file = Pkg.dir("JuliaFEM")*"/test/testdata/primitives.med"
+    mesh = aster_read_mesh(mesh_file, "UNITSQUARE_6X4")
+
+    field = Problem(Heat, "unit square, 6x4 triangular mesh", 1)
+    field.elements = create_elements(mesh, "UNITSQUARE")
+    field.properties.formulation = "2D"
+    update!(field, "temperature thermal conductivity", 1.0)
+    update!(field, "temperature load", -6.0)
+
+    bc = Problem(Dirichlet, "u₀(x,y) = 1 + x² + 2y²", 1, "temperature")
+    #bc.properties.order = 2
+    #bc.properties.dual_basis = true
+    bc.properties.variational = false
+    bc.elements = create_elements(mesh, "FACE1", "FACE2", "FACE3", "FACE4")
+    function u0(element, ip, time)
+        x, y = element("geometry", ip, time)
+        return 1 + x^2 + 2*y^2
+    end
+    update!(bc, "temperature 1", u0)
+
+    solver = LinearSolver(field, bc)
+    solver()
+
+    T_fem = Float64[]
+    T_acc = Float64[]
+    for (nid, X) in field("geometry")
+#       info("$nid -> $X")
+        push!(T_fem, field("temperature", X)[1])
+        push!(T_acc, 1.0 + X[1]^2 + 2*X[2]^2)
+    end
+
+    for element in bc.elements
+        for (X, T_fem) in zip(element("geometry", 0.0), element("temperature", 0.0))
+            x, y = X
+            T_acc = 1.0 + x^2 + 2*y^2
+#           info("(x,y) = ($x,$y), T_acc = $T_acc, T_fem = $T_fem")
+        end
+    end
+
+    @test maximum(abs(T_fem-T_acc)) < 1.0e-12
+end
diff --git a/test/test_postprocess.jl b/test/test_postprocess.jl
new file mode 100644
index 0000000..2d0b6bf
--- /dev/null
+++ b/test/test_postprocess.jl
@@ -0,0 +1,52 @@
+using JuliaFEM
+using JuliaFEM.Test
+
+@testset "get nodal values" begin
+    el1 = Element(Seg2, [1, 2])
+    el2 = Element(Seg2, [2, 3])
+    X = Dict{Int64, Vector{Float64}}(
+        1 => [0.0],
+        2 => [1.0],
+        3 => [2.0])
+    T = Dict{Int64, Vector{Float64}}(
+        1 => [0.0],
+        2 => [1.0],
+        3 => [0.0])
+    P = Problem(Heat, "foo", 1)
+    push!(P, el1, el2)
+    update!(P, "geometry", X)
+    update!(P, "temperature", T)
+    @test isnan(P("temperature", [-0.1]))
+    @test isapprox(P("temperature", [0.0]), [0.0])
+    @test isapprox(P("temperature", [0.5]), [0.5])
+    @test isapprox(P("temperature", [1.0]), [1.0])
+    @test isapprox(P("temperature", [1.5]), [0.5])
+    @test isapprox(P("temperature", [2.0]), [0.0])
+    @test isnan(P("temperature", [ 2.1]))
+end
+
+
+@testset "interpolate from set of elements" begin
+    el1 = Element(Seg2, [1, 2])
+    el2 = Element(Seg2, [2, 3])
+    X = Dict{Int64, Vector{Float64}}(
+        1 => [0.0],
+        2 => [1.0],
+        3 => [2.0])
+    T = Dict{Int64, Vector{Float64}}(
+        1 => [0.0],
+        2 => [1.0],
+        3 => [0.0])
+    P = Problem(Heat, "foo", 1)
+    push!(P, el1, el2)
+    update!(P, "geometry", X)
+    update!(P, "temperature", T)
+    @test isnan(P("temperature", [-0.1]))
+    @test isapprox(P("temperature", [0.0]), [0.0])
+    @test isapprox(P("temperature", [0.5]), [0.5])
+    @test isapprox(P("temperature", [1.0]), [1.0])
+    @test isapprox(P("temperature", [1.5]), [0.5])
+    @test isapprox(P("temperature", [2.0]), [0.0])
+    @test isnan(P("temperature", [ 2.1]))
+end
+
diff --git a/test/test_preprocess_aster_reader.jl b/test/test_preprocess_aster_reader.jl
index 9e23e2d..7276e7c 100644
--- a/test/test_preprocess_aster_reader.jl
+++ b/test/test_preprocess_aster_reader.jl
@@ -151,9 +151,9 @@ end
     @test length(mesh2.elements) == 1
 end
 
-function calculate_volume(eltype::Symbol)
+function calculate_volume(mesh_name::String, eltype::Symbol)
     fn = Pkg.dir("JuliaFEM") * "/test/testdata/primitives.med"
-    mesh = aster_read_mesh(fn, "$eltype")
+    mesh = aster_read_mesh(fn, mesh_name)
     elements = create_elements(mesh, eltype)
     V = 0.0
     time = 0.0
@@ -168,15 +168,22 @@ function calculate_volume(eltype::Symbol)
     return V
 end
 
-@testset "calculate volume for primitives" begin
-    @test isapprox(calculate_volume(:Tet4), 1/6)
-    @test isapprox(calculate_volume(:Tet10), 1/6)
-    @test isapprox(calculate_volume(:Hex8), 2^3)
-    @test isapprox(calculate_volume(:Hex20), 2^3)
-    @test isapprox(calculate_volume(:Hex27), 2^3)
-#   @test isapprox(get_volume("PE6"), V)
-#   @test isapprox(get_volume("PY5"), V)
-#   @test isapprox(get_volume("P15"), V)
-#   @test isapprox(get_volume("P13"), V)
+@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/2))
+#   @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
 
diff --git a/test/testdata/primitives.hdf b/test/testdata/primitives.hdf
index 233bdc1..b7f77fa 100644
Binary files a/test/testdata/primitives.hdf and b/test/testdata/primitives.hdf differ
diff --git a/test/testdata/primitives.med b/test/testdata/primitives.med
index 851c4e3..addaffb 100644
Binary files a/test/testdata/primitives.med and b/test/testdata/primitives.med differ