diff --git a/src/JuliaFEM.jl b/src/JuliaFEM.jl
index 7cca197..22b1cba 100644
--- a/src/JuliaFEM.jl
+++ b/src/JuliaFEM.jl
@@ -115,7 +115,8 @@ include("preprocess_aster_reader.jl")
 export aster_create_elements, parse_aster_med_file, is_aster_mail_keyword,
        parse_aster_header, aster_parse_nodes, aster_renumber_nodes!,
        aster_renumber_elements!, aster_combine_meshes, aster_read_mesh,
-       filter_by_element_set, filter_by_element_id, MEDFile
+       filter_by_element_set, filter_by_element_id, MEDFile, aster_read_data,
+       aster_read_mesh_names, aster_read_node_sets, aster_read_nodes, RMEDFile
 end
 
 function get_mesh(mesh_name::AbstractString, args...; kwargs...)
diff --git a/src/abaqus.jl b/src/abaqus.jl
index a52dfbd..c9817d5 100644
--- a/src/abaqus.jl
+++ b/src/abaqus.jl
@@ -672,7 +672,7 @@ function abaqus_download(name)
         fn = rstrip(ENV["ABAQUS_DOWNLOAD_DIR"], '/') * "/" * fn
     end
     if !isfile(fn)
-        info("Downloading model $name from $url to $fn")
+        info("Downloading model $name ...")
         download("$url/$name.inp", fn)
     end
     return 0
diff --git a/src/postprocess_utils.jl b/src/postprocess_utils.jl
index db37b8f..f6b169b 100644
--- a/src/postprocess_utils.jl
+++ b/src/postprocess_utils.jl
@@ -130,7 +130,7 @@ function calc_nodal_values!(elements::Vector, field_name, field_dim, time;
     end
 end
 
-function calc_nodal_values!(problem::Problem, field_name, field_dim, time)
+function calc_nodal_values!(problem::Problem, field_name::AbstractString, field_dim::Int, time::Float64)
     # after all, it's just a mass matrix ...
 #   isempty(problem.assembly.M) && assemble!(problem, time, Val{:mass_matrix}; density=1.0, dual_basis=false, dim=1)
 #   M = sparse(problem.assembly.M)
@@ -141,7 +141,7 @@ end
 """
 Return node ids + vector of values 
 """
-function get_nodal_vector(elements, field_name, time)
+function get_nodal_vector(elements::Vector, field_name::AbstractString, time::Float64)
     f = Dict()
     for element in elements
         for (c, v) in zip(get_connectivity(element), element[field_name](time))
@@ -337,6 +337,18 @@ function call(problem::Problem, field_name::AbstractString, X::Vector, time::Flo
     return fillna
 end
 
+function call(solver::Solver, field_name::AbstractString, X::Vector, time::Float64; fillna=NaN)
+    for problem in get_problems(solver)
+        for element in get_elements(problem)
+            if inside(element, X, time)
+                xi = get_local_coordinates(element, X, time)
+                return element(field_name, xi, time)
+            end
+        end
+    end
+    return fillna
+end
+
 """ Calculate area of cross-section. """
 function calculate_area(problem::Problem, X=[0.0, 0.0], time=0.0)
     A = 0.0
diff --git a/src/preprocess.jl b/src/preprocess.jl
index 95a0375..e2d9a0d 100644
--- a/src/preprocess.jl
+++ b/src/preprocess.jl
@@ -165,6 +165,18 @@ function reorder_element_connectivity!(mesh::Mesh, mapping::Dict{Symbol, Vector{
     end
 end
 
+function JuliaFEM.Problem{P<:FieldProblem}(mesh::Mesh, ::Type{P}, name::AbstractString, dimension::Int64)
+    problem = Problem{P}(name, dimension, "none", [], Dict(), Assembly(), P())
+    problem.elements = create_elements(mesh, name)
+    return problem
+end
+
+function JuliaFEM.Problem{P<:BoundaryProblem}(mesh::Mesh, ::Type{P}, name, dimension, parent_field_name)
+    problem = Problem{P}(name, dimension, parent_field_name, [], Dict(), Assembly(), P())
+    problem.elements = create_elements(mesh, name)
+    return problem
+end
+
 """
 Swap surface element connectivity s.t. normals point outward
 """
diff --git a/src/preprocess_aster_reader.jl b/src/preprocess_aster_reader.jl
index 2d3e611..e116441 100644
--- a/src/preprocess_aster_reader.jl
+++ b/src/preprocess_aster_reader.jl
@@ -55,19 +55,17 @@ function parse(mesh, ::Type{Val{:CODE_ASTER_MAIL}})
 end
 
 
-"""
-Code Aster binary file (.med), which is exported from SALOME.
-"""
+""" Code Aster binary file (.med). """
 type MEDFile
     data :: Dict
 end
 
 function MEDFile(fn)
-    MEDFile(h5read(fn, "/"))
+    return MEDFile(h5read(fn, "/"))
 end
 
 function get_mesh_names(med::MEDFile)
-    return collect(keys(med.data["FAS"]))
+    return sort(collect(keys(med.data["FAS"])))
 end
 
 function get_nodes(med::MEDFile, nsets, mesh_name)
@@ -91,6 +89,13 @@ end
 
 function get_node_sets(med::MEDFile, mesh_name)
     ns = Dict{Int64, Symbol}(0 => :NALL)
+    if !haskey(med.data["FAS"], mesh_name)
+        warn("Mesh $mesh_name not found from med file.")
+        meshes = get_mesh_names(med)
+        all_meshes = join(meshes, ", ")
+        warn("Available meshes: $all_meshes")
+        error("Mesh $mesh_name not found.")
+    end
     haskey(med.data["FAS"][mesh_name], "NOEUD") || return ns
     nsets = med.data["FAS"][mesh_name]["NOEUD"]
     for nset in keys(nsets)
@@ -250,4 +255,63 @@ function aster_read_mesh(fn, mesh_name=nothing; reorder_element_connectivity=tru
     return mesh
 end
 
-# TODO: refactor and remove obsolete stuff.
+""" Code Aster result file (.rmed). """
+type RMEDFile
+    data :: Dict
+end
+
+function RMEDFile(fn)
+    return RMEDFile(h5read(fn, "/"))
+end
+
+""" Return nodes from result med file. """
+function aster_read_nodes(rmed::RMEDFile)
+    increments = keys(rmed.data["ENS_MAA"]["MAIL"])
+    @assert length(increments) == 1
+    increment = first(increments)
+    nodes = rmed.data["ENS_MAA"]["MAIL"][increment]["NOE"]
+    node_names = nodes["NOM"]
+    node_coords = nodes["COO"]
+    nnodes = length(node_names)
+    dim = round(Int, length(node_coords)/nnodes)
+    node_coords = reshape(node_coords, nnodes, dim)'
+    stripper(node_name) = strip(ascii(pointer(convert(Vector{UInt8}, node_name))))
+    node_names = map(stripper, node_names)
+    # INFO: quite safe assumption is that id is in node name, i.e. N1 => 1, N123 => 123
+    node_id(node_name) = parse(matchall(r"\d+", node_name)[1])
+    node_ids = map(node_id, node_names)
+    nodes = Dict([j => node_coords[:,j] for j in node_ids])
+    return nodes
+end
+
+""" Read nodal field from rmed file. """
+function aster_read_data(rmed::RMEDFile, field_name; field_type=:NODE,
+                         info_fields=true, node_ids=nothing)
+
+    if contains(field_name, "ELGA")
+        field_type = :GAUSS
+    end
+
+    if node_ids == nothing
+        nodes = aster_read_nodes(rmed)
+        node_ids = sort(collect(keys(nodes)))
+    end
+
+    if info_fields
+        field_names = keys(rmed.data["CHA"])
+        all_fields = join(field_names, ", ")
+        info("results: $all_fields")
+    end
+
+    chdata = rmed.data["CHA"]["RESU____$field_name"]
+    @assert length(chdata) == 1
+    increment = chdata[first(keys(chdata))]
+    if field_type == :NODE
+        data = increment["NOE"]["MED_NO_PROFILE_INTERNAL"]["CO"]
+        results = Dict([j => data[j] for j in node_ids])
+    else
+        error("Unable to read result of type $field_type: not implemented")
+    end
+    return results
+end
+
diff --git a/src/problems.jl b/src/problems.jl
index 0e36e7a..9514651 100644
--- a/src/problems.jl
+++ b/src/problems.jl
@@ -104,10 +104,10 @@ julia> prob2 = Problem(Elasticity, 3)
 
 """
 function Problem{P<:FieldProblem}(::Type{P}, name::AbstractString, dimension::Int64)
-    Problem{P}(name, dimension, "none", [], Dict(), Assembly(), P())
+    return Problem{P}(name, dimension, "none", [], Dict(), Assembly(), P())
 end
 function Problem{P<:FieldProblem}(::Type{P}, dimension::Int64)
-    Problem{P}("$P problem", dimension, "none", [], Dict(), Assembly(), P())
+    return Problem{P}("$P problem", dimension, "none", [], Dict(), Assembly(), P())
 end
 
 """ Construct a new boundary problem.
@@ -120,13 +120,13 @@ julia> bc1 = Problem(Dirichlet, "support", 3, "displacement")
 
 """
 function Problem{P<:BoundaryProblem}(::Type{P}, name, dimension, parent_field_name)
-    Problem{P}(name, dimension, parent_field_name, [], Dict(), Assembly(), P())
+    return Problem{P}(name, dimension, parent_field_name, [], Dict(), Assembly(), P())
 end
 function Problem{P<:BoundaryProblem}(::Type{P}, main_problem::Problem)
     name = "$P problem"
     dimension = get_unknown_field_dimension(main_problem)
     parent_field_name = get_unknown_field_name(main_problem)
-    Problem{P}(name, dimension, parent_field_name, [], Dict(), Assembly(), P())
+    return Problem{P}(name, dimension, parent_field_name, [], Dict(), Assembly(), P())
 end
 
 function get_formulation_type{P<:FieldProblem}(problem::Problem{P})
diff --git a/src/problems_mortar_3d.jl b/src/problems_mortar_3d.jl
index 7499ed0..853becf 100644
--- a/src/problems_mortar_3d.jl
+++ b/src/problems_mortar_3d.jl
@@ -86,6 +86,21 @@ function get_cells(P, C)
     info("indices = $indices")
 end
 
+""" Test does P contain q. """
+function contains{T}(P::Vector{T}, q::T; check_is_close=true, rtol=1.0e-5)
+    if q in P
+        return true
+    end
+    if check_is_close
+        for p in P
+            if isapprox(p, q; rtol=rtol)
+                return true
+            end
+        end
+    end
+    return false
+end
+
 function get_polygon_clip(xs, xm, n; debug=false)
     # objective: search does line xm1 - xm2 clip xs
     nm = length(xm)
@@ -103,7 +118,7 @@ function get_polygon_clip(xs, xm, n; debug=false)
     # 2. test is slave point inside master, if yes, add to clip
     for i=1:ns
         if vertex_inside_polygon(xs[i], xm)
-            xs[i] in P && continue
+            contains(P, xs[i]) && continue
             debug && info("2. $(xs[i]) inside M -> push")
             push!(P, xs[i])
         end
@@ -126,7 +141,7 @@ function get_polygon_clip(xs, xm, n; debug=false)
             q = xs1 + t*(xs2 - xs1)
             #info("t=$t, q=$q, q ∈ xm ? $(vertex_inside_polygon(q, xm))")
             if vertex_inside_polygon(q, xm)
-                q in P && continue
+                contains(P, q) && continue
                 debug && info("3. $q inside M -> push")
                 push!(P, q)
             end
diff --git a/src/solvers.jl b/src/solvers.jl
index fe7da6d..ab8431a 100644
--- a/src/solvers.jl
+++ b/src/solvers.jl
@@ -217,6 +217,19 @@ function create_projection(C::SparseMatrixCSC, g; S=nothing, tol=1.0e-12)
     return P, h
 end
 
+""" Assume C is invertible. """
+function create_projection(C, g, ::Type{Val{:invertible}})
+    nz1, nz2 = get_nonzeros(C)
+    P = spzeros(size(C)...)
+    for j=1:size(C,1)
+        j in nz1 && continue
+        P[j,j] = 1.0
+    end
+    v = lufact(C[nz1,nz2]) \ full(g[nz1])
+    return P, v
+end
+
+
 
 """
 Solve linear system using LDLt factorization (SuiteSparse). This version
@@ -332,20 +345,34 @@ function solve_linear_system(solver::Solver; F=nothing, empty_assemblies_before_
 end
 
 """ Default assembler for solver. """
-function assemble!(solver::Solver; show_info=true)
+function assemble!(solver::Solver; show_info=true, timing=true)
     show_info && info("Assembling problems ...")
     t0 = Base.time()
+    assembly_times = Dict()
     nproblems = 0
     ndofs = 0
     for problem in solver.problems
+        t00 = Base.time()
         empty!(problem.assembly)
         assemble!(problem, solver.time)
         nproblems += 1
-        ndofs = max(ndofs, size(problem.assembly.K, 2))
+        Ks = size(problem.assembly.K, 2)
+        Cs = size(problem.assembly.C1, 2)
+        ndofs = max(ndofs, Ks, Cs)
+        t11 = Base.time()
+        assembly_times[problem.name] = t11-t00
     end
     solver.ndofs = ndofs
     t1 = round(Base.time()-t0, 2)
     show_info && info("Assembled $nproblems problems in $t1 seconds. ndofs = $ndofs.")
+    if timing
+        info("Assembly times:")
+        for (i, problem) in enumerate(solver.problems)
+            pn = problem.name
+            pt = round(assembly_times[pn], 2)
+            info("$i $pn $pt")
+        end
+    end
 end
 
 function get_unknown_fields(solver::Solver)
diff --git a/src/solvers_modal.jl b/src/solvers_modal.jl
index a298c2e..6ff8381 100644
--- a/src/solvers_modal.jl
+++ b/src/solvers_modal.jl
@@ -23,7 +23,7 @@ function Modal(nev=10, which=:SM)
     solver = Modal(false, Vector(), Matrix(), nev, which)
 end
 
-function call(solver::Solver{Modal}; show_info=true, debug=false)
+function call(solver::Solver{Modal}; show_info=true, debug=false, bc_invertible=false)
     show_info && info(repeat("-", 80))
     show_info && info("Starting natural frequency solver")
     show_info && info("Increment time t=$(round(solver.time, 3))")
@@ -48,16 +48,59 @@ function call(solver::Solver{Modal}; show_info=true, debug=false)
     if solver.properties.geometric_stiffness
         K += Kg
     end
-    
+
     @assert nnz(D) == 0
     @assert C1 == C2
 
     tic()
-    P, h = create_projection(C1, g)
+
+    if bc_invertible
+        P, h = create_projection(C1, g, Val{:invertible})
+    else
+        P, h = create_projection(C1, g)
+    end
     K_red = P'*K*P
     M_red = P'*M*P
+    # make sure matrices are symmetric
     K_red = 1/2*(K_red + K_red')
     M_red = 1/2*(M_red + M_red')
+
+#=
+    ndim = size(C1,1)
+    nz = get_nonzero_rows(C1)
+    nz = setdiff(collect(1:ndim), nz)
+    g = zeros(ndim)
+    P = spzeros(ndim, ndim)
+    for j in nz
+        P[j,j] = 1.0
+    end
+    K_red = P'*K*P
+    M_red = P'*M*P
+    # make sure matrices are symmetric
+    K_red = 1/2*(K_red + K_red')
+    M_red = 1/2*(M_red + M_red')
+
+    #=
+    K_red = K[nz,nz]
+    M_red = M[nz,nz]
+    # make sure matrices are symmetric
+    K_red = 1/2*(K_red + K_red')
+    M_red = 1/2*(M_red + M_red')
+    =#
+
+    #=
+    K_red = copy(K)
+    M_red = copy(M)
+    for j=1:size(K_red)
+        j in nz && continue
+        K_red[j,:] = 0.0
+        K_red[:,j] = 0.0
+        M_red[j,:] = 0.0
+        M_red[:,j] = 0.0
+    end
+    =#
+=#
+
     t1 = round(toq(), 2)
     info("Eliminated dirichlet boundaries in $t1 seconds.")
 
@@ -79,16 +122,30 @@ function call(solver::Solver{Modal}; show_info=true, debug=false)
         om2, X = eigs(K_red[nz,nz], M_red[nz,nz]; nev=props.nev, which=props.which)
     catch
         info("failed to calculate eigenvalues")
-        info("K sym?", issym(K_red[nz,nz]))
-        info("M sym?", issym(M_red[nz,nz]))
-        info("K posdef?", isposdef(K_red[nz,nz]))
-        info("M posdef?", isposdef(M_red[nz,nz]))
+        info("reduced system")
+        info("is K symmetric? ", issym(K_red[nz,nz]))
+        info("is M symmetric? ", issym(M_red[nz,nz]))
+        info("is K positive definite? ", isposdef(K_red[nz,nz]))
+        info("is M positive definite? ", isposdef(M_red[nz,nz]))
         k1 = maximum(abs(K_red[nz,nz] - K_red[nz,nz]'))
         m1 = maximum(abs(M_red[nz,nz] - M_red[nz,nz]'))
-        info("K skewness ", k1)
-        info("M skewness ", m1)
+        info("K 'skewness' (max(abs(K - K'))) = ", k1)
+        info("M 'skewness' (max(abs(M - M'))) = ", m1)
+
+        info("original matrix")
+        info("is K symmetric? ", issym(K[nz,nz]))
+        info("is M symmetric? ", issym(M[nz,nz]))
+        info("is K positive definite? ", isposdef(K[nz,nz]))
+        info("is M positive definite? ", isposdef(M[nz,nz]))
+        k1 = maximum(abs(K[nz,nz] - K[nz,nz]'))
+        m1 = maximum(abs(M[nz,nz] - M[nz,nz]'))
+        info("K 'skewness' (max(abs(K - K'))) = ", k1)
+        info("M 'skewness' (max(abs(M - M'))) = ", m1)
+
         rethrow()
     end
+    info("Eigenvalues computed in $t1 seconds. Eigenvalues: $om2")
+
     props.eigvals = om2
     props.eigvecs = zeros(ndofs, length(om2))
     v = zeros(ndofs)
@@ -98,7 +155,6 @@ function call(solver::Solver{Modal}; show_info=true, debug=false)
         props.eigvecs[:,i] = P*v + g
     end
     t1 = round(toq(), 2)
-    info("Eigenvalues computed in $t1 seconds. Eigenvalues: $om2")
 
     for i=1:length(om2)
         freq = real(sqrt(om2[i])/(2.0*pi))
diff --git a/src/sparse.jl b/src/sparse.jl
index d382980..fbc6209 100644
--- a/src/sparse.jl
+++ b/src/sparse.jl
@@ -163,6 +163,12 @@ function get_nonzero_columns(A::Union{SparseMatrixCOO, Matrix})
     return get_nonzero_columns(sparse(A))
 end
 
+function get_nonzeros(C::Union{SparseMatrixCSC, Matrix})
+    nz1 = get_nonzero_rows(C)
+    nz2 = get_nonzero_columns(C)
+    return (nz1, nz2)
+end
+
 function size(A::SparseMatrixCOO)
     isempty(A) && return (0, 0)
     return maximum(A.I), maximum(A.J)
@@ -173,7 +179,7 @@ function size(A::SparseMatrixCOO, idx::Int)
 end
 
 """ Matrix norm. Automatically convert to dense when asking for 2-norm for small matrices. """
-function Base.norm(A::SparseMatrixCOO, p=Inf; maxdim=1000)
+function norm(A::SparseMatrixCOO, p=Inf; maxdim=1000)
     dim = size(A, 1)
     if p == 2 && dim > maxdim
         info("Assembly norm: dim = $dim > $maxdim and p=$p, not making dense matrices for operation.")
diff --git a/test/test_heat_3d_two_rings.jl b/test/test_heat_3d_two_rings.jl
new file mode 100644
index 0000000..17ab52e
--- /dev/null
+++ b/test/test_heat_3d_two_rings.jl
@@ -0,0 +1,51 @@
+# 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
+
+#=
+Two rings, RING1 = inner, RING2 = outer, RINGS combined mesh. Set T=1.0 for
+inner ring and T=2.0 for outer ring, measure temperature from middle of ring.
+Results are calculated using Code Aster for comparison.
+=#
+@testset "test 3d heat, two rings, and compare to CA solution" begin
+    meshfile = Pkg.dir("JuliaFEM") * "/test/testdata/primitives.med"
+    mesh = aster_read_mesh(meshfile, "RINGS_UNION")
+
+    rings = Problem(Heat, "RINGS", 1)
+#   rings.elements = create_elements(mesh; element_type=:Tet4)
+    rings.elements = create_elements(mesh, "RING1", "RING2")
+    update!(rings.elements, "temperature thermal conductivity", 1.0)
+    bc_inner = Problem(Dirichlet, "INNER SURFACE", 1, "temperature")
+    bc_inner.elements = create_elements(mesh, "RING1_INNER")
+    bc_outer = Problem(Dirichlet, "OUTER SURFACE", 1, "temperature")
+    bc_outer.elements = create_elements(mesh, "RING2_OUTER")
+    update!(bc_inner, "temperature 1", 1.0)
+    update!(bc_outer, "temperature 1", 2.0)
+    info("# of elements in RING1_INNER = ", length(bc_inner.elements))
+    info("# of elements in RING2_OUTER = ", length(bc_outer.elements))
+    solver = LinearSolver(rings, bc_inner, bc_outer)
+    solver()
+
+    temp_jf = rings("temperature", 0.0)
+
+    fn = Pkg.dir("JuliaFEM") * "/test/testdata/rings.rmed"
+    results = RMEDFile(fn)
+    nodes = aster_read_nodes(results)
+    temp_ca = aster_read_data(results, "TEMP")
+
+    passed = true
+    for j in sort(collect(keys(temp_jf)))
+        X = nodes[j]
+        T1 = temp_jf[j]
+        T2 = temp_ca[j]
+        rtol = norm(T1-T2) / max(T1,T2)
+        @printf "% 5i : %8.5f %8.5f %8.5f | %8.5f %8.5f | %8.5e\n" j X... T1 T2 rtol
+        passed &= rtol < 1.0e-12
+    end
+    @test passed
+end
+
diff --git a/test/test_modal_analysis_elasticity.jl b/test/test_modal_analysis_elasticity.jl
index 6f94fb4..3ca4551 100644
--- a/test/test_modal_analysis_elasticity.jl
+++ b/test/test_modal_analysis_elasticity.jl
@@ -40,10 +40,29 @@ Fixed-fixed solution is ωᵢ = λᵢ²√(EI/ρA) , where λᵢ = cosh(λᵢℓ
 2:  7.853204624095838
 3: 10.995607838001671
 
+Youngs modulus is tuned such that lowest eigenfrequency matches 1.0
+
+5 lowest eigenfrequencies using Code Aster and Tet4 elements:
+numéro    fréquence (HZ)     norme d'erreur
+    1       1.19789E+00        2.20137E-12
+    2       1.20179E+00        1.99034E-12
+    3       3.07391E+00        3.29226E-13
+    4       3.08812E+00        2.91550E-13
+    5       4.87370E+00        2.95986E-13
+
+5 lowest eigenfrequencies using Code Aster and Tet10 elements:
+numéro    fréquence (HZ)     norme d'erreur
+    1       9.65942E-01        1.54950E-11
+    2       9.66160E-01        1.62712E-11
+    3       2.52127E+00        2.06544E-12
+    4       2.52187E+00        1.77970E-12
+    5       3.48584E+00        9.96170E-13
+
+
 [1] De Silva, Clarence W. Vibration: fundamentals and practice. CRC press, 2006, p.355
 =#
 
-@testset "long rod under point load" begin
+@testset "long rod natural frequencies" begin
     mesh_file = Pkg.dir("JuliaFEM") * "/test/testdata/primitives.med"
     mesh = aster_read_mesh(mesh_file, "CYLINDER_20_TET10")
 #   for (id, coords) in mesh.nodes
@@ -51,7 +70,8 @@ Fixed-fixed solution is ωᵢ = λᵢ²√(EI/ρA) , where λᵢ = cosh(λᵢℓ
 #   end
     body = Problem(Elasticity, "rod", 3)
     body.elements = create_elements(mesh, "CYLINDER")
-    E = 50475.44814745859
+    #E = 50475.44814745859
+    E = 50475.5
     rho = 1.0
     update!(body.elements, "youngs modulus", E)
     update!(body.elements, "poissons ratio", 0.3)
@@ -121,13 +141,27 @@ Fixed-fixed solution is ωᵢ = λᵢ²√(EI/ρA) , where λᵢ = cosh(λᵢℓ
     info("freq_a = $freq_a")
 
     solver = Solver(Modal, body, fixed1, fixed2)
+    solver.properties.nev = 5
     solver()
-    freqs = keys(body["displacement"])
+    freqs_jf = sqrt(solver.properties.eigvals)/(2.0*pi)
+    # with Tet4 elements
+    #freqs_ca = [1.19789E+00, 1.20179E+00, 3.07391E+00, 3.08813E+00, 4.87370E+00]
+    # with Tet10 elements
+    freqs_ca = [9.65942E-01, 9.66160E-01, 2.52127E+00, 2.52187E+00, 3.48584E+00]
 
-    rtol1 = norm(freq_sa - freqs[2])/max(freq_sa, freqs[2])
-    rtol2 = norm(freq_a - freqs[2])/max(freq_a, freqs[2])
+    # looks that juliafem results are more close to 1.0, maybe different integration order
+    rtol1 = norm(freq_sa - freqs_jf[1])/max(freq_sa, freqs_jf[1])
+    rtol2 = norm(freq_a - freqs_jf[1])/max(freq_a, freqs_jf[1])
     info("rtol 1 = $rtol1, rtol 2 = $rtol2")
-    @test rtol2 < 1.0e-2
+    passed = true
+    for (f1, f2) in zip(freqs_jf, freqs_ca)
+        rtol = norm(f1-f2) / max(f1,f2)
+        @printf "JF: %8.5e | CA: %8.5e | rtol: %8.5e\n" f1 f2 rtol
+        passed &= (rtol < 3.0e-2)
+    end
+    @test rtol2 < 3.0e-2
+    @test passed
+
 #=
     result = XDMF()
     for (i, freq) in enumerate(freqs)
@@ -138,5 +172,34 @@ Fixed-fixed solution is ωᵢ = λᵢ²√(EI/ρA) , where λᵢ = cosh(λᵢℓ
     end
     xdmf_save!(result, "/tmp/rod_nf.xmf")
 =#
+
+end
+
+@testset "eigenvalues of cube (tet4)" begin
+    meshfile = Pkg.dir("JuliaFEM") * "/test/testdata/primitives.med"
+    mesh = aster_read_mesh(meshfile, "CUBE_TET4")
+    cube = Problem(mesh, Elasticity, "CUBE", 3)
+    update!(cube.elements, "youngs modulus", 10000.0)
+    update!(cube.elements, "poissons ratio", 0.3)
+    update!(cube.elements, "density", 10.0)
+    sym23 = create_elements(mesh, "FACE231")
+    update!(sym23, "displacement 1", 0.0)
+    sym13 = create_elements(mesh, "FACE131")
+    update!(sym13, "displacement 2", 0.0)
+    sym12 = create_elements(mesh, "FACE121")
+    update!(sym12, "displacement 3", 0.0)
+    bcs = Problem(Dirichlet, "bcs", 3, "displacement")
+    bcs.elements = [sym23; sym13; sym12]
+    solver = Solver(Modal)
+    solver.properties.nev = 5
+    push!(solver, cube, bcs)
+    solver()
+    freqs_jf = sqrt(solver.properties.eigvals)/(2.0*pi)
+    freqs_ca = [3.73724E+00, 3.73724E+00, 4.93519E+00, 6.59406E+00, 7.65105E+00]
+    for (f1, f2) in zip(freqs_jf, freqs_ca)
+        rtol = norm(f1-f2) / max(f1,f2)
+        @printf "JF: %8.5e | CA: %8.5e | rtol: %8.5e\n" f1 f2 rtol
+        @test rtol < 1.0e-5
+    end
 end
 
diff --git a/test/test_modal_analysis_elasticity_2.jl b/test/test_modal_analysis_elasticity_2.jl
new file mode 100644
index 0000000..cfded39
--- /dev/null
+++ b/test/test_modal_analysis_elasticity_2.jl
@@ -0,0 +1,68 @@
+# 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 "eigenvalues of CYLINDER1" begin
+    meshfile = Pkg.dir("JuliaFEM") * "/test/testdata/primitives.med"
+    mesh = aster_read_mesh(meshfile, "CYLINDER_1_TET4")
+    cylinder = Problem(mesh, Elasticity, "CYLINDER", 3)
+    update!(cylinder.elements, "youngs modulus", 10000.0)
+    update!(cylinder.elements, "poissons ratio", 0.3)
+    update!(cylinder.elements, "density", 10.0)
+    bc1 = create_elements(mesh, "FACE_YZ1")
+    update!(bc1, "displacement 1", 0.0)
+    update!(bc1, "displacement 2", 0.0)
+    update!(bc1, "displacement 3", 0.0)
+    bcs = Problem(Dirichlet, "bcs", 3, "displacement")
+    bcs.elements = bc1
+    solver = Solver(Modal)
+    solver.properties.nev = 3
+    push!(solver, cylinder, bcs)
+    solver()
+    freqs_jf = sqrt(solver.properties.eigvals)/(2.0*pi)
+    freqs_ca = [4.84532E+00, 4.90698E+00, 8.33813E+00]
+    passed = []
+    for (f1, f2) in zip(freqs_jf, freqs_ca)
+        rtol = norm(f1-f2) / max(f1,f2)
+        @printf "JF: %8.5e | CA: %8.5e | rtol: %8.5e\n" f1 f2 rtol
+        push!(passed, rtol < 1.0e-5)
+    end
+    @test reduce(&, passed)
+end
+
+@testset "eigenvalues of CYLINDER20" begin
+    meshfile = Pkg.dir("JuliaFEM") * "/test/testdata/primitives.med"
+    mesh = aster_read_mesh(meshfile, "CYLINDER_20_TET4")
+    cylinder = Problem(mesh, Elasticity, "CYLINDER", 3)
+    #update!(cylinder.elements, "youngs modulus", 10.0e6)
+    update!(cylinder.elements, "youngs modulus", 50475.5)
+    update!(cylinder.elements, "poissons ratio", 0.3)
+    #update!(cylinder.elements, "density", 10.0)
+    update!(cylinder.elements, "density", 1.0)
+    bc1 = create_elements(mesh, "FACE1", "FACE2")
+    update!(bc1, "displacement 1", 0.0)
+    update!(bc1, "displacement 2", 0.0)
+    update!(bc1, "displacement 3", 0.0)
+    bcs = Problem(Dirichlet, "bcs", 3, "displacement")
+    bcs.elements = bc1
+    solver = Solver(Modal)
+    solver.properties.nev = 3
+    push!(solver, cylinder, bcs)
+    solver()
+    freqs_jf = sqrt(solver.properties.eigvals)/(2.0*pi)
+    #freqs_ca = [8.82848E-01, 8.85353E-01, 5.30286E+00] # only face1 fixed
+    #freqs_ca = [5.33185E+00, 5.34920E+00, 1.36820E+01] # face1 and face2 fixed
+    freqs_ca = [1.19789E+00, 1.20179E+00, 3.07391E+00]
+    passed = []
+    for (f1, f2) in zip(freqs_jf, freqs_ca)
+        rtol = norm(f1-f2) / max(f1,f2)
+        @printf "JF: %8.5e | CA: %8.5e | rtol: %8.5e\n" f1 f2 rtol
+        push!(passed, rtol < 1.0e-5)
+    end
+    @test reduce(&, passed)
+end
+
diff --git a/test/test_mortar_3d_mesh_tie_modal.jl b/test/test_mortar_3d_mesh_tie_modal.jl
new file mode 100644
index 0000000..cfe2804
--- /dev/null
+++ b/test/test_mortar_3d_mesh_tie_modal.jl
@@ -0,0 +1,175 @@
+# 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
+
+#=
+test subjects:
+- modal analysis, with mesh tie contact
+
+Fixed-fixed solution is ωᵢ = λᵢ²√(EI/ρA) , where λᵢ = cosh(λᵢℓ)cos(λᵢℓ)
+
+1:  4.730040744862704
+2:  7.853204624095838
+3: 10.995607838001671
+
+[1] De Silva, Clarence W. Vibration: fundamentals and practice. CRC press, 2006, p.355
+
+Code Aster solution:
+--------------------
+numéro    fréquence (HZ)     norme d'erreur
+    1       1.12946E+00        5.81018E-12
+    2       1.13141E+00        6.33463E-12
+    3       2.93779E+00        6.53408E-13
+    4       2.94143E+00        5.43970E-13
+    5       4.51684E+00        5.43252E-13
+=#
+
+
+comm_CA = """
+DEBUT(PAR_LOT="NON")
+
+MAIL = LIRE_MAILLAGE(FORMAT="MED", NOM_MED="CYLINDER_20_SPLITTED")
+
+MO = AFFE_MODELE(
+    MAILLAGE=MAIL,
+    AFFE=_F(TOUT="OUI",
+            PHENOMENE="MECANIQUE", MODELISATION="3D"))
+
+MAT = DEFI_MATERIAU(
+    ELAS=_F(E=50475.45, NU=0.3, RHO=1.0))
+
+CHMAT = AFFE_MATERIAU(
+    MAILLAGE=MAIL,
+    AFFE=_F(TOUT="OUI", MATER=MAT))
+
+BC1 = AFFE_CHAR_MECA(
+    MODELE=MO,
+    DDL_IMPO=(
+        _F(GROUP_MA=("CYLINDER_20_1_FACE1"), DX=0, DY=0, DZ=0)))
+
+BC2 = AFFE_CHAR_MECA(
+    MODELE=MO,
+    DDL_IMPO=(
+        _F(GROUP_MA=("CYLINDER_20_2_FACE2"), DX=0, DY=0, DZ=0)))
+
+# ESCL = SLAVE
+# MAIT = MASTER
+BC3 = AFFE_CHAR_MECA(
+    MODELE=MO,
+    LIAISON_MAIL=_F(
+        GROUP_MA_ESCL="CYLINDER_20_1_FACE2",
+        GROUP_MA_MAIT="CYLINDER_20_2"))
+
+# assemble material stiffness matrix
+
+RIGEL = CALC_MATR_ELEM(
+    MODELE=MO,
+    OPTION="RIGI_MECA",
+    CHAM_MATER=CHMAT,
+    CHARGE=(BC1, BC2, BC3))
+
+NUMEDDL = NUME_DDL(
+    MATR_RIGI=RIGEL)
+
+RIGAS = ASSE_MATRICE(
+    MATR_ELEM=RIGEL,
+    NUME_DDL=NUMEDDL)
+
+# assemble mass matrix
+
+MASSEL = CALC_MATR_ELEM(
+    MODELE=MO,
+    OPTION="MASS_MECA",
+    CHAM_MATER=CHMAT,
+    CHARGE=(BC1, BC2, BC3))
+
+MASSAS = ASSE_MATRICE(
+    MATR_ELEM=MASSEL,
+    NUME_DDL=NUMEDDL)
+
+# modal analysis, without geometric stiffness
+
+BRESU = CALC_MODES(
+    MATR_RIGI=RIGAS,
+    MATR_MASS=MASSAS,
+    OPTION="BANDE",
+    CALC_FREQ=_F(
+        FREQ=(0.0, 5.0)))
+
+# modal analysis, with geometric stiffness
+
+BRESU = NORM_MODE(
+    reuse=BRESU,
+    MODE=BRESU,
+    NORME="TRAN")
+
+IMPR_RESU(
+    MODELE=MO,
+    FORMAT="RESULTAT",
+    RESU=_F(RESULTAT=BRESU))
+
+IMPR_RESU(
+    FORMAT="MED",
+    UNITE=80,
+    RESU=_F(RESULTAT=BRESU))
+
+FIN()
+"""
+
+@testset "splitted rod with tie contact" begin
+    # CYLINDER_20_1_FACE1 -- CYLINDER_20_1_FACE2 -- CYLINDER_20_2_FACE_1 -- CYLINDER_20_2_FACE_2
+    mesh_file = Pkg.dir("JuliaFEM") * "/test/testdata/primitives.med"
+    mesh = aster_read_mesh(mesh_file, "CYLINDER_20_SPLITTED")
+    body1 = Problem(mesh, Elasticity, "CYLINDER_20_1", 3)
+    body2 = Problem(mesh, Elasticity, "CYLINDER_20_2", 3)
+    for body in [body1, body2]
+        update!(body.elements, "youngs modulus", 54475.45)
+        update!(body.elements, "poissons ratio", 0.3)
+        update!(body.elements, "density", 1.0)
+    end
+    bc1 = Problem(mesh, Dirichlet, "CYLINDER_20_1_FACE1", 3, "displacement")
+    bc2 = Problem(mesh, Dirichlet, "CYLINDER_20_2_FACE2", 3, "displacement")
+    for bc in [bc1, bc2]
+        update!(bc.elements, "displacement 1", 0.0)
+        update!(bc.elements, "displacement 2", 0.0)
+        update!(bc.elements, "displacement 3", 0.0)
+    end
+    interface = Problem(Mortar, "interface between bodies", 3, "displacement")
+    slave = create_elements(mesh, "CYLINDER_20_1_FACE2")
+    master = create_elements(mesh, "CYLINDER_20_2_FACE1")
+    update!(slave, "master elements", master)
+    interface.elements = [slave; master]
+
+    solver = Solver(Modal, body1, body2, bc1, bc2, interface)
+    solver.properties.nev = 5
+    solver.properties.which = :SM
+    solver()
+    freqs_jf = sqrt(solver.properties.eigvals)/(2*pi)
+
+    freqs_ca = [1.12946E+00, 1.13141E+00, 2.93779E+00, 2.94143E+00, 4.51684E+00]
+    freq_jf = freqs_jf[1]
+    freq_ca = freqs_ca[1]
+    rtol = norm(freq_jf - freq_ca)/max(freq_jf, freq_ca)
+    info("rtol = $rtol")
+    for (i, freq) in enumerate(freqs_jf)
+        @printf "mode %i | freq JuliaFEM %8.3f | freq Code Aster %8.3f\n" i freqs_jf[i] freqs_ca[i]
+    end
+    if rtol > 1.0e-3
+        outfile = tempname() * ".xmf"
+        info("Something went wrong, results are saved to $outfile")
+        result = XDMF()
+        elems = [body1.elements; body2.elements]
+        for (i, freq) in enumerate(freqs_jf)
+            xdmf_new_result!(result, elems, freq)
+            xdmf_save_field!(result, elems, freq, "displacement"; field_type="Vector")
+        end
+        xdmf_save!(result, outfile)
+    end
+    @test rtol < 0.05
+
+end
+
diff --git a/test/test_mortar_3d_mesh_tie_two_rings.jl b/test/test_mortar_3d_mesh_tie_two_rings.jl
new file mode 100644
index 0000000..d6c1172
--- /dev/null
+++ b/test/test_mortar_3d_mesh_tie_two_rings.jl
@@ -0,0 +1,61 @@
+# 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
+
+#=
+Two rings, RING1 is inner, RING2 is outer. Inner diameter is from 0.8 .. 0.9 and
+outer ring is 0.9 .. 1.0. Contact surface pair is RING1_OUTER <- RING2_INNER.
+Put constant temperature 1.0 for inner surface of inner ring and 2.0 for outer
+surface of outer ring. We should expect constant temperature in contact surface.
+This is conforming mesh so result should match to the conforming situation.
+=#
+@testset "test that curved interface transfers constant field without error, two rings problem" begin
+    meshfile = Pkg.dir("JuliaFEM") * "/test/testdata/primitives.med"
+    mesh = aster_read_mesh(meshfile, "RINGS")
+
+    ring1 = Problem(Heat, "RING1", 1)
+    ring1.elements = create_elements(mesh, "RING1")
+    update!(ring1.elements, "temperature thermal conductivity", 1.0)
+
+    ring2 = Problem(Heat, "RING2", 1)
+    ring2.elements = create_elements(mesh, "RING2")
+    update!(ring2.elements, "temperature thermal conductivity", 1.0)
+
+    bc_inner = Problem(Dirichlet, "INNER SURFACE", 1, "temperature")
+    bc_inner.elements = create_elements(mesh, "RING1_INNER")
+    update!(bc_inner, "temperature 1", 1.0)
+
+    bc_outer = Problem(Dirichlet, "OUTER SURFACE", 1, "temperature")
+    bc_outer.elements = create_elements(mesh, "RING2_OUTER")
+    update!(bc_outer, "temperature 1", 2.0)
+
+    interface = Problem(Mortar, "interface between rings", 1, "temperature")
+    interface_slave = create_elements(mesh, "RING1_OUTER")
+    interface_master = create_elements(mesh, "RING2_INNER")
+    interface.elements = [interface_slave; interface_master]
+    update!(interface_slave, "master elements", interface_master)
+
+    solver = LinearSolver(ring1, ring2, bc_inner, bc_outer, interface)
+    solver()
+
+    fn = Pkg.dir("JuliaFEM") * "/test/testdata/rings.rmed"
+    results = RMEDFile(fn)
+    nodes = aster_read_nodes(results)
+    temp_ca = aster_read_data(results, "TEMP")
+
+    passed = true
+    for j in sort(collect(keys(nodes)))
+        X = nodes[j]
+        T1 = solver("temperature", X, 0.0)
+        T2 = temp_ca[j]
+        rtol = norm(T1-T2) / max(T1,T2)
+        @printf "% 5i : %8.5f %8.5f %8.5f | %8.5f %8.5f | %8.5f\n" j X... T1 T2 rtol
+        passed = passed && (rtol < 1.0e-12)
+    end
+    @test passed
+end
+
diff --git a/test/test_preprocess_aster_reader.jl b/test/test_preprocess_aster_reader.jl
index d6f72b1..ba01619 100644
--- a/test/test_preprocess_aster_reader.jl
+++ b/test/test_preprocess_aster_reader.jl
@@ -147,3 +147,12 @@ end
 #   @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_projection.jl b/test/test_projection.jl
index 2902e43..d7af575 100644
--- a/test/test_projection.jl
+++ b/test/test_projection.jl
@@ -3,6 +3,7 @@
 
 using JuliaFEM
 using JuliaFEM.Testing
+using JuliaFEM.Preprocess
 
 @testset "test projection" begin
     C = [
@@ -28,5 +29,66 @@ using JuliaFEM.Testing
     h_expected = [1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
     @test isapprox(full(P), P_expected)
     @test isapprox(full(h), h_expected)
+    nz = [5, 6, 7, 8]
+    info("is P'P positive definite? ", isposdef(P[:,nz]'P[:,nz]))
+    @test isposdef(P[:,nz]'P[:,nz])
+end
+
+@testset "test creating projection matrix from invertible problem" begin
+    # simple 3 element poisson problem
+    k = [1.0 -1.0; -1.0 1.0]
+    K = zeros(4, 4)
+    K[1:2,1:2] += k
+    K[2:3,2:3] += k
+    K[3:4,3:4] += k
+    # first dof homogeneous bc, last dof u₄ = 1
+    C = zeros(4, 4)
+    C[1,1] = 1.0
+    C[4,4] = 1.0
+    g = zeros(4)
+    g[1] = 0.0
+    g[4] = 1.0
+    P, h = create_projection(sparse(C), g, Val{:invertible})
+    info("P = ")
+    dump(full(P))
+    P_expected = zeros(4, 4)
+    P_expected[2,2] = P_expected[3,3] = 1.0
+    @test isapprox(full(P), P_expected)
+    #h_expected = [0.5, 1.0]
+    #@test isapprox(h, h_expected)
+end
+
+@testset "projection between surfaces" begin
+    # FIXME: creating mortar projection takes very long time.
+    meshfile = Pkg.dir("JuliaFEM") * "/test/testdata/joint.med"
+    isfile(meshfile) || return
+    mesh = aster_read_mesh(meshfile, "JOINT")
+
+    # top
+    bc1 = Problem(Dirichlet, "fixed1", 3, "displacement")
+    bc1.elements = create_elements(mesh, "FIXED1")
+    update!(bc1.elements, "displacement 1", 0.0)
+    update!(bc1.elements, "displacement 2", 0.0)
+    update!(bc1.elements, "displacement 3", 0.0)
+
+    # bottom
+    bc2 = Problem(Dirichlet, "fixed2", 3, "displacement")
+    bc2.elements = create_elements(mesh, "FIXED2")
+    update!(bc2.elements, "displacement 1", 0.0)
+    update!(bc2.elements, "displacement 2", 0.0)
+    update!(bc2.elements, "displacement 3", 0.0)
+
+    # joint
+    contact = Problem(Mortar, "joint", 3, "displacement")
+    master_elements = create_elements(mesh, "BODY1_TO_BODY2")
+    slave_elements = create_elements(mesh, "BODY2_TO_BODY1")
+    update!(slave_elements, "master elements", master_elements)
+    contact.elements = [master_elements; slave_elements]
+
+    solver = Solver(Modal)
+    push!(solver, bc1, bc2, contact)
+    assemble!(solver)
+    Kb, C1, C2, D, fb, g = get_boundary_assembly(solver)
+    P, h = create_projection(C1, g)
 end
 
diff --git a/test/testdata/joint.med b/test/testdata/joint.med
new file mode 100644
index 0000000..6d79496
Binary files /dev/null and b/test/testdata/joint.med differ
diff --git a/test/testdata/primitives.hdf b/test/testdata/primitives.hdf
index 8d1ec86..d27c7ce 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 7f6195a..8393351 100644
Binary files a/test/testdata/primitives.med and b/test/testdata/primitives.med differ
diff --git a/test/testdata/rings.rmed b/test/testdata/rings.rmed
new file mode 100644
index 0000000..e0f1a0d
Binary files /dev/null and b/test/testdata/rings.rmed differ