diff --git a/src/JuliaFEM.jl b/src/JuliaFEM.jl
index a797fef..b4a4812 100644
--- a/src/JuliaFEM.jl
+++ b/src/JuliaFEM.jl
@@ -1,14 +1,112 @@
 # This file is a part of JuliaFEM.
 # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
 
-# __precompile__()
-
 """
-This is JuliaFEM -- Finite Element Package
+    JuliaFEM.jl - an open source solver for both industrial and academia usage
+
+The JuliaFEM software library is a framework that allows for the distributed
+processing of large Finite Element Models across clusters of computers using
+simple programming models. It is designed to scale up from single servers to
+thousands of machines, each offering local computation and storage. The basic
+design principle is: everything is nonlinear. All physics models are nonlinear
+from which the linearization are made as a special cases.
+
+# Examples
+
+Typical workflow to use JuliaFEM to solve a partial differential equations,
+is 1) read mesh, 2) create elements, 3) update element properties, 4) create
+problems, 5) create analysis, and 6) run analysis. A simple linear static
+analysis of elastic block is clarifying these steps.
+
+```julia
+using JuliaFEM
+```
+
+1. Usually the first thing to do is to create a geometry of domain. Typically
+   this is done by reading a mesh file from disk. Currently JuliaFEM supports
+   reading a mesh from Code Aster file format (using `aster_read_mesh`) and
+   from ABAQUS file format (using `abaqus_read_mesh`).
+
+```julia
+mesh = aster_read_mesh("mesh.med")
+```
+
+2. Next step is to create one or several sets of elements from a mesh. This is
+   done using a function `create_elements(mesh, set_name)`.
+
+```julia
+body_elements = create_elements(mesh, "body")
+traction_elements = create_elements(mesh, "traction")
+bc_elements = create_elements(mesh, "bc")
+```
+
+3. In JuliaFEM, all properties of elements are given using so called fields.
+   Fields can depend from time or spatial coordinates elements. In special cases
+   field is constant in time, spatial or both directions. Updating fields to
+   element is done using `update!`-function.
+
+```julia
+update!(body_elements, "youngs modulus", 210.0e3)
+update!(body_elements, "poissons ratio", 0.3)
+update!(traction_elements, "surface pressure", 100.0)
+update!(bc_elements, "displacement 1", 0.0)
+update!(bc_elements, "displacement 2", 0.0)
+update!(bc_elements, "displacement 3", 0.0)
+```
+
+4. The physics considered to be solve is given using `Problem`. Problem type can
+   be e.g. `Elasticity` for hyperelasticity, `Heat` for solving the heat equation
+   and so on. Problems are defined by giving the problem type as first argument,
+   problem name in second argument and the last argument is giving the dimension
+   of problem, meaning degrees of freedom connected to each node. After problems
+   are created, elements are added to them by using function `add_elements!`.
+
+```julia
+body = Problem(Elasticity, "body", 3)
+traction = Problem(Elasticity, "traction", 3)
+bc = Problem(Dirichlet, "bc", 3, "displacement")
+add_elements!(body, body_elements)
+add_elements!(traction, traction_elements)
+add_elements!(bc, bc_elements)
+```
+
+5. After geometry and physics is defined, we next define what kind of analysis
+   are we going to perform. Analysis can be, for example, quasistatic analysis,
+   analysis of dynamics of system, natural frequency analysis and so on. Some
+   other special analysis types also exists, like performing model dimension
+   reduction by creating super-elements or running optimization loop, given
+   geometry, another analysis and initial conditions. For simplicity, we now
+   create a linear quasistatic analysis of given problems. Problems are added
+   to analysis using `add_problems!`.
+
+```julia
+analysis = Analysis(Linear)
+add_problems!(analysis, body, traction, bc)
+```
+
+6. The last thing to do is to request the results of analysis to be written to
+   disk for later use and actually perform the analysis. Currently, Xdmf output
+   is supported, which can then be read using ParaView.
+
+```julia
+xdmf_output = Xdmf("analysis_results")
+add_results_writer!(analysis, xdmf_output)
+run!(analysis)
+close(xdmf)
+```
+
+After analysis is ready, types and variables can be accessed using REPL or
+Jupyter notebook for further postprocessing. Simulation can also be written
+into a function to be a part of a larger analysis process. For more information
+about JuliaFEM, please visit our website at
+
+    www.juliafem.org
+
 """
 module JuliaFEM
 
-using Reexport
+using SparseArrays, LinearAlgebra, Statistics
+using Reexport, ForwardDiff, LightXML, HDF5
 
 @reexport using FEMBase
 import FEMBase: get_unknown_field_name, get_unknown_field_dimension,
@@ -21,12 +119,7 @@ export @timeit, print_timer
 import Base: getindex, setindex!, convert, length, size, isapprox,
              similar, start, first, next, done, last, endof, vec,
              ==, +, -, *, /, haskey, copy, push!, isempty, empty!,
-             append!, sparse, full, read
-
-module Testing
-using Base.Test
-export @test, @testset, @test_throws
-end
+             append!, read, copy
 
 using AbaqusReader
 using AsterReader
@@ -70,31 +163,29 @@ include("problems_contact_3d.jl")
 #include("problems_contact_3d_autodiff.jl")
 export Contact
 
-# Preprocess module
-
 module Preprocess
-using FEMBase
+end
+
+using FEMBase, SparseArrays, LinearAlgebra
 include("preprocess.jl")
-export create_elements, Mesh, add_node!, add_nodes!, add_element!,
-       add_elements!, add_element_to_element_set!, add_node_to_node_set!,
+export create_elements, Mesh, add_node!, add_nodes!,
+       add_element_to_element_set!, add_node_to_node_set!,
        find_nearest_nodes, find_nearest_node, reorder_element_connectivity!,
        create_node_set_from_element_set!, filter_by_element_set
 include("preprocess_abaqus_reader.jl")
 export abaqus_read_mesh, create_surface_elements, create_nodal_elements
 include("preprocess_aster_reader.jl")
 export aster_read_mesh
-end
 
 # Postprocess module
 
 module Postprocess
-using FEMBase
-using FEMBase: get_elements
+end
+
 include("postprocess_utils.jl")
 export calc_nodal_values!, get_nodal_vector, get_nodal_dict, copy_field!,
-       calculate_area, calculate_center_of_mass,
-       calculate_second_moment_of_mass, extract
-end
+       calculate_area, calculate_center_of_mass, calculate_second_moment_of_mass,
+       extract
 
 include("deprecations.jl")
 
diff --git a/src/deprecations.jl b/src/deprecations.jl
index 29de2b4..92998c8 100644
--- a/src/deprecations.jl
+++ b/src/deprecations.jl
@@ -10,5 +10,5 @@ function assemble!(problem::Problem, element::Element, time=0.0)
 end
 
 module Abaqus
-using JuliaFEM.Preprocess: create_surface_elements
+using JuliaFEM: create_surface_elements
 end
diff --git a/src/io.jl b/src/io.jl
index 6b93e0e..1316a41 100644
--- a/src/io.jl
+++ b/src/io.jl
@@ -36,7 +36,7 @@ function Xdmf(name::String; version="3.0", overwrite=false)
 
     if isfile(h5file)
         if overwrite
-            info("Result file $h5file exists, removing old file.")
+            @debug("Result file $h5file exists, removing old file.")
             rm(h5file)
         else
             error("Result file $h5file exists, use Xdmf($name; overwrite=true) to rewrite results")
@@ -45,7 +45,7 @@ function Xdmf(name::String; version="3.0", overwrite=false)
 
     if isfile(xmlfile)
         if overwrite
-            info("Result file $xmlfile exists, removing old file.")
+            @debug("Result file $xmlfile exists, removing old file.")
             rm(xmlfile)
         else
             error("Result file $xmlfile exists, use Xdmf($name; overwrite=true) to rewrite results")
@@ -209,9 +209,9 @@ function traverse(xdmf::Xdmf, x::XMLElement, attr_name::String)
         items = split(attr_name, '/')
         new_item = xdmf_filter(childs, first(items))
         if new_item == nothing
-            info("traverse: childs:")
+            @debug("traverse: childs:")
             for child in childs
-                info(LightXML.name(child))
+                @debug(LightXML.name(child))
             end
             error("traverse: failed, items = $items, xdmf_filter not find child")
         end
@@ -261,9 +261,13 @@ function save!(xdmf::Xdmf)
     save_file(doc, xmffile(xdmf))
 end
 
+function Base.close(xdmf::Xdmf)
+    close(xdmf.hdf)
+end
+
 function new_dataitem(xdmf::Xdmf, path::String, data::Array{T,N}) where {T,N}
     dataitem = new_element("DataItem")
-    datatype = replace("$T", "64", "")
+    datatype = replace("$T", "64" => "")
     dimensions = join(reverse(size(data)), " ")
     set_attribute(dataitem, "DataType", datatype)
     set_attribute(dataitem, "Dimensions", dimensions)
@@ -271,7 +275,7 @@ function new_dataitem(xdmf::Xdmf, path::String, data::Array{T,N}) where {T,N}
     if xdmf.format == "HDF"
         hdf = basename(h5file(xdmf))
         if exists(xdmf.hdf, path)
-            info("Xdmf: $path already existing in h5 file, not overwriting.")
+            @debug("Xdmf: $path already existing in h5 file, not overwriting.")
         else
             write(xdmf.hdf, path, data)
         end
@@ -324,6 +328,28 @@ global const xdmf_element_mapping = Dict(
         "Wedge15" => "Wedge_15",
         "Hex20" => "Hex_20")
 
+"""
+    get_spatial_collection()
+
+Return a SpatialCollection at given time either by creating new one or returning
+existing one.
+"""
+function get_spatial_collection(temporal_collection, time)
+    for spatial_collection in get_elements_by_tagname(temporal_collection, "Grid")
+        time_element = find_element(spatial_collection, "Time")
+        time_value = Meta.parse(attribute(time_element, "Value"; required=true))
+        isapprox(time_value, time) && return spatial_collection
+    end
+    # did not find, create new one
+    spatial_collection = new_child(temporal_collection, "Grid")
+    set_attribute(spatial_collection, "GridType", "Collection")
+    set_attribute(spatial_collection, "Name", "Problems")
+    set_attribute(spatial_collection, "CollectionType", "Spatial")
+    time_element = new_child(spatial_collection, "Time")
+    set_attribute(time_element, "Value", time)
+    return spatial_collection
+end
+
 """
     update_xdmf!(xdmf, problem, time, fields)
 
@@ -337,20 +363,20 @@ julia> update_xdmf!(p1, 0.0, ["displacement", "temperature"])
 """
 function update_xdmf!(xdmf::Xdmf, problem::Problem, time::Float64, fields::Vector)
 
-    info("Xdmf: storing fields $fields of problem $(problem.name) at time $time")
+    @debug("Xdmf: storing fields $fields of problem $(problem.name) at time $time")
 
     # 1. find domain
     xml = xdmf.xml
     domain = find_element(xml, "Domain")
     if domain == nothing
-        info("Xdmf: Domain not found, creating.")
+        @debug("Xdmf: Domain not found, creating.")
         domain = new_child(xml, "Domain")
     end
 
     # 2. find for TemporalCollection
     temporal_collection = find_element(domain, "Grid")
     if temporal_collection == nothing
-        info("Xdmf: Temporal collection not found, creating.")
+        @debug("Xdmf: Temporal collection not found, creating.")
         temporal_collection = new_child(domain, "Grid")
         set_attribute(temporal_collection, "GridType", "Collection")
         set_attribute(temporal_collection, "Name", "Time")
@@ -361,43 +387,18 @@ function update_xdmf!(xdmf::Xdmf, problem::Problem, time::Float64, fields::Vecto
     collection_type = attribute(temporal_collection, "CollectionType"; required=true)
     @assert collection_type == "Temporal"
 
-    # 3. find for SpatialCollection at given time
-    spatial_collection = nothing
-    spatial_collection_exists = false
-    for spatial_collection in get_elements_by_tagname(temporal_collection, "Grid")
-        time_element = find_element(spatial_collection, "Time")
-        time_value = parse(attribute(time_element, "Value"; required=true))
-        if isapprox(time_value, time)
-            info("Xdmf: SpatialCollection for time $time already exists.")
-            spatial_collection_exists = true
-            break
-        end
-    end
-
-    if !spatial_collection_exists
-        info("Xdmf: SpatialCollection for time $time not found, creating.")
-        spatial_collection = new_child(temporal_collection, "Grid")
-        set_attribute(spatial_collection, "GridType", "Collection")
-        set_attribute(spatial_collection, "Name", "Problems")
-        set_attribute(spatial_collection, "CollectionType", "Spatial")
-        time_element = new_child(spatial_collection, "Time")
-        set_attribute(time_element, "Value", time)
-    end
-
-    # 3.1 make sure that Grid element we found really is SpatialCollection
-    collection_type = attribute(spatial_collection, "CollectionType"; required=true)
-    @assert collection_type == "Spatial"
+    spatial_collection = get_spatial_collection(temporal_collection, time)
 
     for frame in get_elements_by_tagname(spatial_collection, "Grid")
         frame_name = attribute(frame, "Name")
         if frame_name == problem.name
-            warn("Xdmf: Already found Grid with name $frame_name for time $time, skipping.")
+            @warn("Xdmf: Already found Grid with name $frame_name for time $time, skipping.")
             return
         end
     end
 
     frame_name = problem.name
-    info("Xdmf: Creating Grid for problem $frame_name")
+    @debug("Xdmf: Creating Grid for problem $frame_name")
     frame = new_child(spatial_collection, "Grid")
     set_attribute(frame, "Name", frame_name)
 
@@ -408,7 +409,7 @@ function update_xdmf!(xdmf::Xdmf, problem::Problem, time::Float64, fields::Vecto
     X_array = hcat([X_dict[nid] for nid in node_ids]...)
     ndim, nnodes = size(X_array)
     geom_type = (ndim == 2 ? "XY" : "XYZ")
-    info("Xdmf: Creating geometry, type = $geom_type, number of nodes = $nnodes")
+    @debug("Xdmf: Creating geometry, type = $geom_type, number of nodes = $nnodes")
     X_dataitem = new_dataitem(xdmf, X_array)
     geometry = new_child(frame, "Geometry")
     set_attribute(geometry, "Type", geom_type)
@@ -419,12 +420,12 @@ function update_xdmf!(xdmf::Xdmf, problem::Problem, time::Float64, fields::Vecto
     nelements = length(all_elements)
     element_types = unique(map(get_element_type, all_elements))
     nelement_types = length(element_types)
-    info("Xdmf: Saving topology of $nelements elements total, $nelement_types different element types.")
+    @debug("Xdmf: Saving topology of $nelements elements total, $nelement_types different element types.")
 
     for element_type in element_types
         elements = collect(filter_by_element_type(element_type, all_elements))
         nelements = length(elements)
-        info("Xdmf: $nelements elements of type $element_type")
+        @debug("Xdmf: $nelements elements of type $element_type")
         sort!(elements, by=get_element_id)
         element_ids = map(get_element_id, elements)
         element_conn = map(element -> [node_mapping[j]-1 for j in get_connectivity(element)], elements)
@@ -446,24 +447,24 @@ function update_xdmf!(xdmf::Xdmf, problem::Problem, time::Float64, fields::Vecto
         @assert node_ids == field_node_ids
         field_dim = length(field_dict[first(field_node_ids)])
         if field_dim == 2
-            info("Xdmf: Field dimension = 2, extending to 3")
+            @debug("Xdmf: Field dimension = 2, extending to 3")
             for nid in field_node_ids
                 field_dict[nid] = [field_dict[nid]; 0.0]
             end
             field_dim = 3
         end
         field_type = Dict(1 => "Scalar", 3 => "Vector", 6 => "Tensor6")[field_dim]
-        info("Xdmf: Saving field $field_name, type = $field_type, dimension = $field_dim, center = $field_center")
+        @debug("Xdmf: Saving field $field_name, type = $field_type, dimension = $field_dim, center = $field_center")
 
         field_array = hcat([field_dict[nid] for nid in field_node_ids]...)
         field_dataitem = new_dataitem(xdmf, field_array)
         attribute = new_child(frame, "Attribute")
-        set_attribute(attribute, "Name", ucfirst(field_name))
+        set_attribute(attribute, "Name", uppercasefirst(field_name))
         set_attribute(attribute, "Center", field_center)
         set_attribute(attribute, "AttributeType", field_type)
         add_child(attribute, field_dataitem)
     end
 
     save!(xdmf)
-    info("Xdmf: all done.")
+    @debug("Xdmf: all done.")
 end
diff --git a/src/materials_plasticity.jl b/src/materials_plasticity.jl
index 9ab0379..500fa0f 100644
--- a/src/materials_plasticity.jl
+++ b/src/materials_plasticity.jl
@@ -27,7 +27,7 @@ function equivalent_stress(stress, ::Type{Val{:type_3d}})
     stress_ten = [stress[1] stress[6] stress[5];
                   stress[6] stress[2] stress[4];
                   stress[5] stress[4] stress[3]]
-    stress_dev = stress_ten - 1/3 * trace(stress_ten) * eye(3)
+    stress_dev = stress_ten - 1/3 * tr(stress_ten) * eye(3)
     s = vec(stress_dev)
     return sqrt(3/2 * dot(s, s))
 end
diff --git a/src/postprocess_utils.jl b/src/postprocess_utils.jl
index 971185c..cf21683 100644
--- a/src/postprocess_utils.jl
+++ b/src/postprocess_utils.jl
@@ -1,12 +1,6 @@
 # This file is a part of JuliaFEM.
 # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
 
-using JuliaFEM
-#using DataFrames
-using HDF5
-using LightXML
-using Formatting
-
 """
 Calculate field values to nodal points from Gauss points using least-squares fitting.
 """
@@ -27,7 +21,7 @@ function calc_nodal_values!(elements::Vector, field_name, field_dim, time;
         A = sparse(A)
         nz = get_nonzero_rows(A)
         A = 1/2*(A + A')
-        F = ldltfact(A[nz,nz])
+        F = ldlt(A[nz,nz])
     end
 
     if b == nothing
@@ -36,7 +30,7 @@ function calc_nodal_values!(elements::Vector, field_name, field_dim, time;
             gdofs = get_connectivity(element)
             for ip in get_integration_points(element)
                 if !haskey(ip, field_name)
-                    info("warning: integration point does not have field $field_name")
+                    @warn("integration point does not have field $field_name")
                     continue
                 end
                 detJ = element(ip, time, Val{:detJ})
@@ -81,8 +75,30 @@ function get_nodal_vector(elements::Vector, field_name::AbstractString, time::Fl
     return node_ids, field
 end
 
-""" Interpolate field from a set of elements. """
-function (problem::Problem)(field_name::AbstractString, X::Vector, time::Float64=0.0; fillna=NaN)
+"""
+    problem(field_name, X, time)
+
+Interpolate field from a set of elements defined in problem. Here, `X` is the
+location inside domain described by elements.
+
+Internally, function loops through all the elements, finding the one containing
+the point `X`. After that, using inverse isoparametric mapping, first find
+dimensionless coordinates (ξ,η,ζ) of that element corresponding to the location
+of point `X` and after that interpolate the values of field under investigation.
+Algorithm can be expected to be somewhat slow for big models, but for tests
+models the performance is good.
+
+# Examples
+
+Having a problem called `body`, one can query the field `displacement` at
+position `X = (1.0, 2.0, 3.0)` and time `t = 1.0`, with the command
+```julia
+X = (1.0, 2.0, 3.0)
+time = 1.0
+u = body("displacement", X, time)
+```
+"""
+function (problem::Problem)(field_name, X, time; fillna=NaN)
     for element in get_elements(problem)
         if inside(element, X, time)
             xi = get_local_coordinates(element, X, time)
@@ -92,8 +108,7 @@ function (problem::Problem)(field_name::AbstractString, X::Vector, time::Float64
     return fillna
 end
 
-""" Interpolate field from a set of elements. """
-function (problem::Problem)(field_name::AbstractString, X::Vector, time::Float64, ::Type{Val{:Grad}}; fillna=NaN)
+function (problem::Problem)(field_name, X, time, ::Type{Val{:Grad}}; fillna=NaN)
     for element in get_elements(problem)
         if inside(element, X, time)
             xi = get_local_coordinates(element, X, time)
@@ -134,7 +149,7 @@ https://en.wikipedia.org/wiki/Center_of_mass
 """
 function calculate_center_of_mass(problem::Problem, X=[0.0, 0.0, 0.0], time=0.0)
     M = 0.0
-    Xc = zeros(X)
+    Xc = zero(X)
     for element in get_elements(problem)
         for ip in get_integration_points(element)
             w = ip.weight*element(ip, time, Val{:detJ})
diff --git a/src/preprocess.jl b/src/preprocess.jl
index d759b9c..2698986 100644
--- a/src/preprocess.jl
+++ b/src/preprocess.jl
@@ -12,10 +12,6 @@
 - etc only topology related stuff
 =#
 
-import Base: copy
-
-using JuliaFEM
-
 mutable struct Mesh
     nodes :: Dict{Int, Vector{Float64}}
     node_sets :: Dict{Symbol, Set{Int}}
@@ -42,8 +38,9 @@ function Mesh(m::Dict)
     mesh.nodes = m["nodes"]
     mesh.elements = m["elements"]
     mesh.element_types = m["element_types"]
-    mesh.surface_sets = m["surface_sets"]
-    mesh.surface_types = m["surface_types"]
+    for (k, v) in m["surface_types"]
+        mesh.surface_types[Symbol(k)] = v
+    end
     for (nset_name, node_ids) in m["node_sets"]
         mesh.node_sets[Symbol(nset_name)] = Set(node_ids)
     end
@@ -100,7 +97,7 @@ the set names to be inserted in the function.
 function create_node_set_from_element_set!(mesh::Mesh, set_names::String...)
     for set_name in set_names
         set_name = Symbol(set_name)
-        info("Creating node set $set_name from element set")
+        @info("Creating node set $set_name from element set")
         node_ids = Set{Int}()
         for elid in mesh.element_sets[set_name]
             push!(node_ids, mesh.elements[elid]...)
@@ -125,9 +122,10 @@ end
 Add an element into the mesh. ´elid´ is the element id, ´eltype´ is the type of
 the element and ´connectivity´ is the connectivity of the element.
 """
-function add_element!(mesh::Mesh, elid::Int, eltype::Symbol, connectivity::Vector{Int})
+function FEMBase.add_element!(mesh::Mesh, elid, eltype, connectivity)
     mesh.elements[elid] = connectivity
     mesh.element_types[elid] = eltype
+    return nothing
 end
 
 """
@@ -135,10 +133,11 @@ end
 
 Add elements into the mesh.
 """
-function add_elements!(mesh::Mesh, elements::Dict{Int, Tuple{Symbol, Vector{Int}}})
+function FEMBase.add_elements!(mesh::Mesh, elements::Dict{Int, Tuple{Symbol, Vector{Int}}})
     for (elid, (eltype, elcon)) in elements
         add_element!(mesh, elid, eltype, elcon)
     end
+    return nothing
 end
 
 """
@@ -157,9 +156,9 @@ end
 """
     copy(mesh)
 
-Copy the mesh.
+Return a copy of the mesh.
 """
-function copy(mesh::Mesh)
+function Base.copy(mesh::Mesh)
     mesh2 = Mesh()
     mesh2.nodes = copy(mesh.nodes)
     mesh2.node_sets = copy(mesh.node_sets)
@@ -208,11 +207,6 @@ function create_element(mesh::Mesh, id::Int)
     return element
 end
 
-"""
-    create_elements(mesh, element_type=nothing)
-
-Create elements from the mesh filtered by their type.
-"""
 function create_elements(mesh::Mesh; element_type=nothing)
     element_ids = collect(keys(mesh.elements))
     if element_type != nothing
@@ -237,12 +231,34 @@ function create_elements(mesh::Mesh, element_sets::Symbol...; element_type=nothi
     end
 
     elements = [create_element(mesh, id) for id in element_ids]
+
+    nelements = length(elements)
+    content = Dict{Symbol, Int}()
+    for elid in element_ids
+        eltype = mesh.element_types[elid]
+        content[eltype] = get(content, eltype, 0) + 1
+    end
+    s = join(("$v x $k" for (k, v) in content), ", ")
+    v = join(element_sets, ", ")
+    @info("Created $nelements elements ($s) from element set: $v.")
+
     return elements
 end
 
-function create_elements(mesh::Mesh, element_sets::AbstractString...; element_type=nothing)
-    element_sets = map(parse, element_sets)
-    return create_elements(mesh, element_sets...; element_type=element_type)
+"""
+    create_elements(mesh::Mesh, element_set::String)
+
+# Examples
+
+Suppose that there is a `mesh` with element set `Body_1`. Creating elements
+based on that element set is done then
+
+```julia
+create_elements(mesh, "Body_1")
+```
+"""
+function create_elements(mesh::Mesh, element_sets::String...)
+    return create_elements(mesh, map(Symbol, element_sets)...)
 end
 
 
diff --git a/src/preprocess_aster_reader.jl b/src/preprocess_aster_reader.jl
index e2fb61e..b5a7444 100644
--- a/src/preprocess_aster_reader.jl
+++ b/src/preprocess_aster_reader.jl
@@ -48,11 +48,12 @@ const med_element_names = Dict{Symbol, Symbol}(
     :P13 => :Pyr13)
 
 """
-    aster_read_mesh(filename::String, mesh_name=nothing; reorder_element_connectivity=true)
+    aster_read_mesh(filename, mesh_name=nothing; reorder_element_connectivity=true)
 
-Read code aster mesh from file and return Mesh instance. If mesh file contains
-several meshes, a name of mesh must be given. By default elements are reordered
-so that they match to the conventions used in JuliaFEM.
+Read code aster mesh from file and return `Mesh` structure.
+
+If mesh file contains several meshes, a name of mesh must be given. By default,
+elements are reordered so that they match to the conventions used in JuliaFEM.
 """
 function aster_read_mesh(filename::String, mesh_name=nothing; reorder_element_connectivity=true)
     m = AsterReader.aster_read_mesh(filename, mesh_name)
@@ -63,5 +64,19 @@ function aster_read_mesh(filename::String, mesh_name=nothing; reorder_element_co
     if reorder_element_connectivity
         reorder_element_connectivity!(mesh, med_connectivity)
     end
+    nnodes = length(mesh.nodes)
+    nelements = length(mesh.elements)
+    @info("Mesh parsed from Code Aster file $filename.")
+    @info("Mesh contains $nnodes nodes and $nelements elements.")
+    for (elset_name, elset_elids) in mesh.element_sets
+        content = Dict{Symbol, Int}()
+        for elid in elset_elids
+            eltype = mesh.element_types[elid]
+            content[eltype] = get(content, eltype, 0) + 1
+        end
+        s = join(("$v x $k" for (k, v) in content), ", ")
+        nels = length(elset_elids)
+        @info("Element set $elset_name contains $nels elements ($s).")
+    end
     return mesh
 end
diff --git a/src/problems_contact_3d.jl b/src/problems_contact_3d.jl
index cf485de..0761fe6 100644
--- a/src/problems_contact_3d.jl
+++ b/src/problems_contact_3d.jl
@@ -4,8 +4,8 @@
 const ContactElements3D = Union{Tri3,Tri6,Quad4,Quad8,Quad9}
 
 function create_orthogonal_basis(n)
-    I = eye(3)
-    k = indmax([norm(cross(n,I[:,k])) for k in 1:3])
+    I = [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0]
+    k = argmax([norm(cross(n,I[:,k])) for k in 1:3])
     t1 = cross(n, I[:,k])/norm(cross(n, I[:,k]))
     t2 = cross(n, t1)
     return t1, t2
@@ -106,7 +106,6 @@ function create_contact_segmentation(slave_element, master_elements, x0, n0, tim
     return result
 end
 
-"Assemble linear surface element to contact problem. """
 function assemble!(problem::Problem{Contact}, slave_element::Element{Tri3}, time::Float64)
 
     props = problem.properties
@@ -134,7 +133,7 @@ function assemble!(problem::Problem{Contact}, slave_element::Element{Tri3}, time
         return
     end
 
-    Ae = eye(nsl)
+    Ae = Matrix{Float64}(I, nsl, nsl)
 
     if problem.properties.dual_basis # construct dual basis
 
@@ -154,7 +153,7 @@ function assemble!(problem::Problem{Contact}, slave_element::Element{Tri3}, time
                     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))
+                    De += w*Matrix(Diagonal(vec(N1)))
                     Me += w*N1'*N1
                 end # integration points done
 
@@ -163,7 +162,7 @@ function assemble!(problem::Problem{Contact}, slave_element::Element{Tri3}, time
         end # master elements done
 
         Ae = De*inv(Me)
-    
+
     end
 
     # loop all polygons
@@ -193,18 +192,18 @@ function assemble!(problem::Problem{Contact}, slave_element::Element{Tri3}, 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 = interpolate(N1, map(+,X1,u1))
                 x_m = interpolate(N2, map(+,X2,u2))
                 ge += w*vec((x_m-x_s)*Phi')
-            
+
             end # integration points done
 
         end # integration cells done
@@ -220,7 +219,7 @@ function assemble!(problem::Problem{Contact}, slave_element::Element{Tri3}, time
             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)
@@ -250,16 +249,16 @@ function assemble!(problem::Problem{Contact}, slave_element::Element{Tri6}, time
                 alp 0.0 alp 0.0 0.0 1.0-2*alp
             ]
     else
-        T = eye(6)
+        T = Matrix(1.0*I, 6, 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)
+    Ae = Matrix(1.0*I, nsl, nsl)
 
     if problem.properties.dual_basis # construct dual basis
 
@@ -329,7 +328,7 @@ function assemble!(problem::Problem{Contact}, slave_element::Element{Tri6}, time
                             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)
+                            De += w*Matrix(Diagonal(N1))
                             Me += w*N1*N1'
                         end # integration points done
 
@@ -342,7 +341,7 @@ function assemble!(problem::Problem{Contact}, slave_element::Element{Tri6}, time
         end # sub slave elements done
 
         Ae = De*inv(Me)
-    
+
     end
 
     # split slave element to linear sub-elements and loop
@@ -352,13 +351,13 @@ function assemble!(problem::Problem{Contact}, slave_element::Element{Tri6}, time
         nsl = length(sub_slave_element)
         X1 = sub_slave_element("geometry", time)
         n1 = sub_slave_element("normal", time)
-            
+
         # create auxiliary plane
         xi = get_mean_xi(sub_slave_element)
         N = vec(get_basis(sub_slave_element, xi, time))
         x0 = interpolate(N, X1)
         n0 = interpolate(N, n1)
-            
+
         # project slave nodes to auxiliary plane
         S = Vector[project_vertex_to_auxiliary_plane(p, x0, n0) for p in X1]
 
@@ -416,7 +415,7 @@ function assemble!(problem::Problem{Contact}, slave_element::Element{Tri6}, 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))
@@ -424,13 +423,13 @@ function assemble!(problem::Problem{Contact}, slave_element::Element{Tri6}, time
 
                         De += w*Phi*N1'
                         Me += w*Phi*N2'
-               
+
                         us = slave_element("displacement", time)
                         um = master_element("displacement", time)
                         xs = interpolate(N1, map(+,Xs,us))
                         xm = interpolate(N2, map(+,Xs,um))
                         ge += w*vec((xm-xs)*Phi')
-            
+
                     end # integration points done
 
                 end # integration cells done
@@ -446,7 +445,7 @@ function assemble!(problem::Problem{Contact}, slave_element::Element{Tri6}, time
                     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)
@@ -513,8 +512,8 @@ function assemble!(problem::Problem{Contact}, time::Float64, ::Type{Val{2}}, ::T
     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)
+    g = Vector(problem.assembly.g, ndofs)
+    c = Vector(problem.assembly.c, ndofs)
 
     maxdim = maximum(size(C1))
     if problem.properties.alpha != 0.0
@@ -550,13 +549,13 @@ function assemble!(problem::Problem{Contact}, time::Float64, ::Type{Val{2}}, ::T
         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)
+        d[get_nonzero_rows(T)] .= 0.0
+        T += sparse(Diagonal(d))
+        invT += sparse(Diagonal(d))
         #invT2 = sparse(inv(full(T)))
-        #info("invT == invT2? ", invT == invT2)
+        #@info("invT == invT2? ", invT == invT2)
         #maxabsdiff = maximum(abs(invT - invT2))
-        #info("max diff = $maxabsdiff")
+        #@info("max diff = $maxabsdiff")
         C1 = C1*invT
         C2 = C2*invT
     end
@@ -572,7 +571,7 @@ function assemble!(problem::Problem{Contact}, time::Float64, ::Type{Val{2}}, ::T
 
     state = problem.properties.contact_state_in_first_iteration
     if problem.properties.iteration == 1
-        info("First contact iteration, initial contact state = $state")
+        @info("First contact iteration, initial contact state = $state")
 
         if state == :AUTO
             avg_gap = mean([weighted_gap[j][1] for j in S])
@@ -582,7 +581,7 @@ function assemble!(problem::Problem{Contact}, time::Float64, ::Type{Val{2}}, ::T
             else
                 state = :UNKNOWN
             end
-            info("Average weighted gap = $avg_gap, std gap = $std_gap, automatically determined contact state = $state")
+            @info("Average weighted gap = $avg_gap, std gap = $std_gap, automatically determined contact state = $state")
         end
 
     end
@@ -602,7 +601,7 @@ function assemble!(problem::Problem{Contact}, time::Float64, ::Type{Val{2}}, ::T
             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.0
             is_inactive[j] = 0
             is_active[j] = 1
@@ -624,7 +623,7 @@ function assemble!(problem::Problem{Contact}, time::Float64, ::Type{Val{2}}, ::T
             is_stick[j] = 0
         end
     end
-    
+
     if (problem.properties.iteration == 1) && (state == :INACTIVE)
         for j in S
             is_inactive[j] = 1
@@ -634,35 +633,31 @@ function assemble!(problem::Problem{Contact}, time::Float64, ::Type{Val{2}}, ::T
         end
     end
 
-    info("# | active | stick | slip | gap | pres | comp")
+    @info("# | active | stick | slip | gap | pres | comp")
     for j in S
         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)
+        str2 = "$(round(weighted_gap[j][1]; digits=3)) | $(round(contact_pressure[j][1]; digits=3)) | $(round(complementarity_condition[j][1]; digits=3))"
+        @info(str1 * str2)
     end
-    
-    # remove inactive nodes from assembly
+
+
     for j in S
         dofs = [3*(j-1)+1, 3*(j-1)+2, 3*(j-1)+3]
+        tdofs = [3*(j-1)+2, 3*(j-1)+3]
         if is_inactive[j] == 1
-            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)
-            C2[tdofs,:] = 0.0
-            g[tdofs] = 0.0
+            # remove inactive nodes from assembly
+            C1[dofs,:] .= 0.0
+            C2[dofs,:] .= 0.0
+            D[dofs,:] .= 0.0
+            g[dofs,:] .= 0.0
+        elseif (is_active[j] == 1) && (is_slip[j] == 1)
+            # constitutive modelling in tangent direction, frictionless contact
+            C2[tdofs,:] .= 0.0
+            g[tdofs] .= 0.0
             normal = normals[j]
             tangent1, tangent2 = create_orthogonal_basis(normal)
-            D[tdofs[1], dofs] = tangent1
-            D[tdofs[2], dofs] = tangent2
+            D[tdofs[1], dofs] .= tangent1
+            D[tdofs[2], dofs] .= tangent2
         end
     end
 
diff --git a/src/problems_dirichlet.jl b/src/problems_dirichlet.jl
index ecab62a..441e7d9 100644
--- a/src/problems_dirichlet.jl
+++ b/src/problems_dirichlet.jl
@@ -24,7 +24,7 @@ function get_dualbasis(element::Element, time::Float64, order=1)
         detJ = element(ip, time, Val{:detJ})
         w = ip.weight*detJ
         N = element(ip, time)
-        De += w*diagm(vec(N))
+        De += w*Matrix(Diagonal(vec(N)))
         Me += w*N'*N
     end
     return De, Me, De*inv(Me)
@@ -38,23 +38,23 @@ function assemble!(problem::Problem{Dirichlet}, time::Float64=0.0;
                    auto_initialize=true)
     # FIXME: boilerplate
     if !isempty(problem.assembly)
-        warn("Assemble problem $(problem.name): problem.assembly is not empty and assembling, are you sure you know what are you doing?")
+        @warn("Assemble problem $(problem.name): problem.assembly is not empty and assembling, are you sure you know what are you doing?")
     end
     if isempty(problem.elements)
-        warn("Assemble problem $(problem.name): problem.elements is empty, no elements in problem?")
+        @warn("Assemble problem $(problem.name): problem.elements is empty, no elements in problem?")
     else
         first_element = first(problem.elements)
         unknown_field_name = get_unknown_field_name(problem)
         if !haskey(first_element, unknown_field_name)
-            warn("Assemble problem $(problem.name): seems that problem is uninitialized.")
+            @warn("Assemble problem $(problem.name): seems that problem is uninitialized.")
             if auto_initialize
-                info("Initializing problem $(problem.name) at time $time automatically.")
+                @info("Initializing problem $(problem.name) at time $time automatically.")
                 initialize!(problem, time)
             end
         end
     end
 
-    if method_exists(assemble_prehook!, Tuple{typeof(problem), Float64})
+    if hasmethod(assemble_prehook!, Tuple{typeof(problem), Float64})
         assemble_prehook!(problem, time)
     end
 
@@ -88,13 +88,13 @@ function assemble!(problem::Problem{Dirichlet}, time::Float64=0.0;
             end
         end
         for (k, v) in field_vals
-            push!(problem.assembly.C1, k, k, 1.0)
-            push!(problem.assembly.C2, k, k, 1.0)
-            push!(problem.assembly.g, k, 1, v)
+            add!(problem.assembly.C1, k, k, 1.0)
+            add!(problem.assembly.C2, k, k, 1.0)
+            add!(problem.assembly.g, k, 1, v)
         end
     end
 
-    if method_exists(assemble_posthook!, Tuple{typeof(problem), Float64})
+    if hasmethod(assemble_posthook!, Tuple{typeof(problem), Float64})
         assemble_posthook!(problem, time)
     end
 end
@@ -112,7 +112,7 @@ function assemble!(assembly::Assembly, problem::Problem{Dirichlet},
     if problem.properties.dual_basis
         De, Me, Ae = get_dualbasis(element, time)
     else
-        Ae = eye(nnodes)
+        Ae = I
         De = zeros(nnodes, nnodes)
         for ip in get_integration_points(element, props.order)
             N = element(ip, time)
diff --git a/src/problems_elasticity.jl b/src/problems_elasticity.jl
index 1f8931a..faff70b 100644
--- a/src/problems_elasticity.jl
+++ b/src/problems_elasticity.jl
@@ -261,9 +261,9 @@ function assemble!(assembly::Assembly,
             :plastic_strain in props.store_fields && update!(ip, "plastic_strain", time => plastic_strain)
 
             #Km += w*BL'*Dtan*BL
-            At_mul_B!(Bt_mul_D, BL, Dtan)
-            A_mul_B!(Bt_mul_D_mul_B, Bt_mul_D, BL)
-            scale!(Bt_mul_D_mul_B, w)
+            mul!(Bt_mul_D, transpose(BL), Dtan)
+            mul!(Bt_mul_D_mul_B, Bt_mul_D, BL)
+            rmul!(Bt_mul_D_mul_B, w)
             for i=1:ndofs^2
                 @inbounds Km[i] += Bt_mul_D_mul_B[i]
             end
@@ -301,8 +301,8 @@ function assemble!(assembly::Assembly,
             end
 
             # internal load
-            At_mul_B!(Bt_mul_S, BL, stress_vec)
-            scale!(Bt_mul_S, w)
+            mul!(Bt_mul_S, transpose(BL), stress_vec)
+            rmul!(Bt_mul_S, w)
             for i=1:ndofs
                 @inbounds f_int[i] += Bt_mul_S[i]
             end
@@ -407,7 +407,7 @@ function get_stress_tensor(problem, element, ip, time)
     nu = element("poissons ratio", ip, time)
     mu = E/(2.0*(1.0+nu))
     la = E*nu/((1.0+nu)*(1.0-2.0*nu))
-    S = la*trace(eps)*I + 2.0*mu*eps
+    S = la*tr(eps)*I + 2.0*mu*eps
     return S
 end
 
@@ -445,13 +445,13 @@ function lsq_fit(problem, elements, field, time)
     A = sparse(A)
     b = sparse(b)
     A = 1/2*(A + A')
-    
+
     nz = get_nonzero_rows(A)
-    F = ldltfact(A[nz,nz])
+    F = ldlt(A[nz,nz])
 
     x = F \ b[nz, :]
 
-    nodal_values = Dict(node_id => vec(full(x[idx,:])) for (idx, node_id) in enumerate(nz))
+    nodal_values = Dict(node_id => Vector(x[idx, :]) for (idx, node_id) in enumerate(nz))
     return nodal_values
 end
 
diff --git a/src/problems_elasticity_2d.jl b/src/problems_elasticity_2d.jl
index 3e80076..78921d6 100644
--- a/src/problems_elasticity_2d.jl
+++ b/src/problems_elasticity_2d.jl
@@ -56,7 +56,7 @@ function assemble(problem::Problem{Elasticity},
 
         if props.finite_strain
             strain = 1/2*(gradu + gradu' + gradu'*gradu)
-            F = eye(dim) + gradu
+            F = I + gradu
             for i=1:size(dN, 2)
                 BL[1, 2*(i-1)+1] += F[1,1]*dN[1,i]
                 BL[1, 2*(i-1)+2] += F[2,1]*dN[1,i]
@@ -67,7 +67,7 @@ function assemble(problem::Problem{Elasticity},
             end
         else # linearized strain
             strain = 1/2*(gradu + gradu')
-            F = eye(dim)
+            F = I
             for i=1:size(dN, 2)
                 BL[1, 2*(i-1)+1] = dN[1,i]
                 BL[2, 2*(i-1)+2] = dN[2,i]
diff --git a/src/problems_mortar.jl b/src/problems_mortar.jl
index 6ec769e..733d6b5 100644
--- a/src/problems_mortar.jl
+++ b/src/problems_mortar.jl
@@ -54,12 +54,12 @@ end
 
 function assemble!(problem::Problem{Mortar}, time::Float64)
     if length(problem.elements) == 0
-        warn("No elements defined in interface $(problem.name), this will result empty assembly!")
+        @warn("No elements defined in interface $(problem.name), this will result empty assembly!")
         return
     end
     if problem.properties.dimension == -1
         problem.properties.dimension = dim = size(first(problem.elements), 1)
-        info("Assuming dimension of mesh tie surface is $dim. If this is wrong set is manually using problem.properties.dimension")
+        @info("Assuming dimension of mesh tie surface is $dim. If this is wrong set is manually using problem.properties.dimension")
     end
     dimension = Val{problem.properties.dimension}
     use_forwarddiff = Val{problem.properties.use_forwarddiff}
@@ -99,7 +99,7 @@ end
 
 """ Function to print useful debug information from interface to find bugs. """
 function diagnose_interface(problem::Problem{Mortar}, time::Float64)
-    info("Diagnosing Mortar interface...")
+    @info("Diagnosing Mortar interface...")
     props = problem.properties
     field_dim = get_unknown_field_dimension(problem)
     field_name = get_parent_field_name(problem)
@@ -108,13 +108,13 @@ function diagnose_interface(problem::Problem{Mortar}, time::Float64)
     I_area = 0.0
 
     if props.split_quadratic_slave_elements
-        info("props.split_quadratic_slave_elements = true")
+        @info("props.split_quadratic_slave_elements = true")
         if !props.linear_surface_elements
-            warn("Mortar3D: split_quadratic_surfaces = true and linear_surface_elements = false maybe have unexpected behavior")
+            @warn("Mortar3D: split_quadratic_surfaces = true and linear_surface_elements = false maybe have unexpected behavior")
         end
         slave_elements = split_quadratic_elements(slave_elements, time)
     end
-    info("Number of slave elements in interface: $(length(slave_elements))")
+    @info("Number of slave elements in interface: $(length(slave_elements))")
 
     # 1. calculate nodal normals and tangents for slave element nodes j ∈ S
     normals = calculate_normals(slave_elements, time, Val{2};
@@ -127,25 +127,25 @@ function diagnose_interface(problem::Problem{Mortar}, time::Float64)
 
     for slave_element in slave_elements
 
-        info(repeat("-", 80))
-        info("Processing slave element $(slave_element.id), type = $(get_element_type(slave_element))")
-        info(repeat("-", 80))
+        @info(repeat("-", 80))
+        @info("Processing slave element $(slave_element.id), type = $(get_element_type(slave_element))")
+        @info(repeat("-", 80))
 
         S_area = 0.0
         S_area_in_contact = 0.0
         for ip in get_integration_points(slave_element)
             S_area += ip.weight*slave_element(ip, time, Val{:detJ})
         end
-        info("Total area of slave element = $S_area")
+        @info("Total area of slave element = $S_area")
 
 
         if props.linear_surface_elements
-            info("Converting slave element to linear surface element")
+            @info("Converting slave element to linear surface element")
             slave_element = convert_to_linear_element(slave_element)
         end
 
         slave_element_nodes = get_connectivity(slave_element)
-        info("Slave element connectivity = $slave_element_nodes")
+        @info("Slave element connectivity = $slave_element_nodes")
         nsl = length(slave_element)
         X1 = slave_element("geometry", time)
         n1 = tuple(collect(normals[j] for j in slave_element_nodes)...)
@@ -155,10 +155,10 @@ function diagnose_interface(problem::Problem{Mortar}, time::Float64)
         N = vec(get_basis(slave_element, xi, time))
         x0 = interpolate(N,X1)
         n0 = interpolate(N,n1)
-        info("Auxiliary plane x0 = $x0, n0 = $n0")
+        @info("Auxiliary plane x0 = $x0, n0 = $n0")
         S = Vector[project_vertex_to_auxiliary_plane(X1[i], x0, n0) for i=1:nsl]
         check_orientation!(S, n0)
-        info("Slave element $(slave_element.id) vertices in auxiliary plane: $S")
+        @info("Slave element $(slave_element.id) vertices in auxiliary plane: $S")
 
         # 3. loop all master elements
         master_elements = slave_element("master elements", time)
@@ -190,17 +190,17 @@ function diagnose_interface(problem::Problem{Mortar}, time::Float64)
                     continue
                 end
                 if length(P) == 1
-                    info("length(P) == 1, shared vertex")
+                    @info("length(P) == 1, shared vertex")
                 end
                 if length(P) == 2
-                    info("length(P) == 2, shared edge")
+                    @info("length(P) == 2, shared edge")
                 end
                 continue
             end
-            info("Master element $(master_element.id) vertices in auxiliary plane = $M")
+            @info("Master element $(master_element.id) vertices in auxiliary plane = $M")
             check_orientation!(P, n0)
             P_area_ = calculate_polygon_area(P)
-            info("Polygon clip found, P=$P, N_P = $(length(P)), area of polygon = $P_area_")
+            @info("Polygon clip found, P=$P, N_P = $(length(P)), area of polygon = $P_area_")
             if isapprox(P_area_, 0.0)
                 error("Polygon P has zero area: $P_area_")
             end
@@ -208,11 +208,11 @@ function diagnose_interface(problem::Problem{Mortar}, time::Float64)
             P_area = 0.0
 
             C0 = calculate_centroid(P)
-            info("Centroid of polygon = $C0")
+            @info("Centroid of polygon = $C0")
 
             # 4. loop integration cells
             all_cells = get_cells(P, C0)
-            info("Polygon is splitted to $(length(all_cells)) integration cells.")
+            @info("Polygon is splitted to $(length(all_cells)) integration cells.")
             for (cell_id, cell) in enumerate(all_cells)
                 C_area = 0.0
                 virtual_element = Element(Tri3, Int[])
@@ -229,7 +229,7 @@ function diagnose_interface(problem::Problem{Mortar}, time::Float64)
                     xi_m, alpha = project_vertex_to_surface(x_gauss, x0, n0, master_element, X2, time)
                     C_area += w
                 end # integration points done
-                info("Cell $cell_id has area of $C_area")
+                @info("Cell $cell_id has area of $C_area")
                 P_area += C_area
                 push!(C_areas, C_area)
             end # integration cells done
@@ -245,15 +245,15 @@ function diagnose_interface(problem::Problem{Mortar}, time::Float64)
 
         S_perc = S_area_in_contact / S_area * 100.0
         push!(S_areas, S_area_in_contact)
-        info("Area of slave element in contact: $S_area_in_contact, it's $S_perc % of total element area")
+        @info("Area of slave element in contact: $S_area_in_contact, it's $S_perc % of total element area")
 
         I_area += S_area_in_contact
 
     end # slave elements done, contact virtual work ready
 
-    info("Area of interface: $I_area")
-    info("Smallest cell area: $(minimum(C_areas))")
-    info("Smallest polygon area: $(minimum(P_areas))")
-    info("Smallest slave element area in contact: $(minimum(S_areas))")
+    @info("Area of interface: $I_area")
+    @info("Smallest cell area: $(minimum(C_areas))")
+    @info("Smallest polygon area: $(minimum(P_areas))")
+    @info("Smallest slave element area in contact: $(minimum(S_areas))")
 
 end
diff --git a/src/problems_mortar_3d.jl b/src/problems_mortar_3d.jl
index d89d9a8..be59949 100644
--- a/src/problems_mortar_3d.jl
+++ b/src/problems_mortar_3d.jl
@@ -37,10 +37,10 @@ function vertex_inside_polygon(q, P; atol=1.0e-3)
         #try
         angle += acos(cosa)
         #catch
-        #    info("Unable to calculate acos($(ForwardDiff.get_value(cosa))) when determining is a vertex inside polygon.")
-        #    info("Polygon is: $(ForwardDiff.get_value(P)) and vertex under consideration is $(ForwardDiff.get_value(q))")
-        #    info("Polygon corner point in loop: A=$(ForwardDiff.get_value(A)), B=$(ForwardDiff.get_value(B))")
-        #    info("c = ||A||*||B|| = $(ForwardDiff.get_value(c))")
+        #    @info("Unable to calculate acos($(ForwardDiff.get_value(cosa))) when determining is a vertex inside polygon.")
+        #    @info("Polygon is: $(ForwardDiff.get_value(P)) and vertex under consideration is $(ForwardDiff.get_value(q))")
+        #    @info("Polygon corner point in loop: A=$(ForwardDiff.get_value(A)), B=$(ForwardDiff.get_value(B))")
+        #    @info("c = ||A||*||B|| = $(ForwardDiff.get_value(c))")
         #    rethrow()
         #end
     end
@@ -128,18 +128,18 @@ function get_polygon_clip(xs::Vector{T}, xm::Vector{T}, n::T) where T
         # 2. find possible intersection
         xm1 = xm[i]
         xm2 = xm[mod(i,nm)+1]
-        #info("intersecting line $xm1 -> $xm2")
+        # @info("intersecting line $xm1 -> $xm2")
         for j=1:ns
             xs1 = xs[j]
             xs2 = xs[mod(j,ns)+1]
-            #info("clipping polygon edge $xs1 -> $xs2")
+            # @info("clipping polygon edge $xs1 -> $xs2")
             tnom = dot(cross(xm1-xs1, xm2-xm1), n)
             tdenom = dot(cross(xs2-xs1, xm2-xm1), n)
             isapprox(tdenom, 0) && continue
             t = tnom/tdenom
             (0 <= t <= 1) || continue
             q = xs1 + t*(xs2 - xs1)
-            #info("t=$t, q=$q, q ∈ xm ? $(vertex_inside_polygon(q, xm))")
+            # @info("t=$t, q=$q, q ∈ xm ? $(vertex_inside_polygon(q, xm))")
             if vertex_inside_polygon(q, xm)
                 approx_in(q, P) && continue
                 push!(P, q)
@@ -180,27 +180,27 @@ function project_vertex_to_surface(p, x0, n0,
         end
     end
     #=
-    info("failed to project vertex from auxiliary plane back to surface")
-    info("element type: $E")
-    info("element connectivity: $(get_connectivity(element))")
-    info("auxiliary plane: x0 = $x0, n0 = $n0")
-    info("element geometry: $(x.data)")
-    info("vertex to project: $p")
-    info("parameter vector before giving up: $theta")
-    info("increment in parameter vector before giving up: $dtheta")
-    info("norm(dtheta) before giving up: $(norm(dtheta))")
-    info("f([0.0, 0.0, 0.0]) = $(f([0.0, 0.0, 0.0]))")
-    info("L([0.0, 0.0, 0.0]) = $(L([0.0, 0.0, 0.0]))")
+    @info("failed to project vertex from auxiliary plane back to surface")
+    @info("element type: $E")
+    @info("element connectivity: $(get_connectivity(element))")
+    @info("auxiliary plane: x0 = $x0, n0 = $n0")
+    @info("element geometry: $(x.data)")
+    @info("vertex to project: $p")
+    @info("parameter vector before giving up: $theta")
+    @info("increment in parameter vector before giving up: $dtheta")
+    @info("norm(dtheta) before giving up: $(norm(dtheta))")
+    @info("f([0.0, 0.0, 0.0]) = $(f([0.0, 0.0, 0.0]))")
+    @info("L([0.0, 0.0, 0.0]) = $(L([0.0, 0.0, 0.0]))")
 
-    info("iterations:")
+    @info("iterations:")
     theta = zeros(3)
     dtheta = zeros(3)
     for i=1:max_iterations
-        info("iter $i, theta = $theta")
-        info("f = $(f(theta))")
-        info("L = $(L(theta))")
+        @info("iter $i, theta = $theta")
+        @info("f = $(f(theta))")
+        @info("L = $(L(theta))")
         dtheta = L(theta) * f(theta)
-        info("dtheta = $(dtheta)")
+        @info("dtheta = $(dtheta)")
         theta -= dtheta
     end
     =#
@@ -273,8 +273,8 @@ function check_orientation!(P, n)
     sort!(P, lt=(A, B) -> begin
         A_proj = Q'*(A-C)
         B_proj = Q'*(B-C)
-        a = atan2(A_proj[3], A_proj[2])
-        b = atan2(B_proj[3], B_proj[2])
+        a = atan(A_proj[3], A_proj[2])
+        b = atan(B_proj[3], B_proj[2])
         return a > b
     end)
 end
@@ -330,7 +330,7 @@ function split_quadratic_elements(elements::Vector, time::Float64)
     n1 = length(elements)
     n2 = length(new_elements)
     if n1 != n2
-        info("Splitted $n1 elements to $n2 (linear) sub-elements")
+        @info("Splitted $n1 elements to $n2 (linear) sub-elements")
     end
     return new_elements
 end
@@ -361,7 +361,7 @@ References
 [Popp2013] Popp, Alexander, et al. "Improved robustness and consistency of 3D contact algorithms based on a dual mortar approach." Computer Methods in Applied Mechanics and Engineering 264 (2013): 67-80.
 """
 function assemble!(problem::Problem{Mortar}, slave_element::Element{E}, time::Real; first_slave_element=false) where E<:Union{Tri3, Quad4}
-    
+
     props = problem.properties
     field_dim = get_unknown_field_dimension(problem)
     field_name = get_parent_field_name(problem)
@@ -385,7 +385,7 @@ function assemble!(problem::Problem{Mortar}, slave_element::Element{E}, time::Re
 
         De = zeros(nsl, nsl)
         Me = zeros(nsl, nsl)
- 
+
         for master_element in master_elements
 
             master_element_nodes = get_connectivity(master_element)
@@ -405,7 +405,7 @@ function assemble!(problem::Problem{Mortar}, slave_element::Element{E}, time::Re
             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)
-                info("Polygon P has zero area: $P_area")
+                @info("Polygon P has zero area: $P_area")
                 continue
             end
 
@@ -421,19 +421,19 @@ function assemble!(problem::Problem{Mortar}, slave_element::Element{E}, time::Re
                     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))
+                    De += w*Matrix(Diagonal(vec(N1)))
                     Me += w*N1'*N1
                 end
             end # integration cells done
 
         end # master elements done
-                    
+
         Ae = De*inv(Me)
 
-        info("Dual basis coefficient matrix: $Ae")
+        @info("Dual basis coefficient matrix: $Ae")
 
     else
-        Ae = eye(nsl)
+        Ae = Matrix(1.0I, nsl, nsl)
     end
 
     for master_element in master_elements
@@ -455,7 +455,7 @@ function assemble!(problem::Problem{Mortar}, slave_element::Element{E}, time::Re
         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)
-            info("Polygon P has zero area: $P_area")
+            @info("Polygon P has zero area: $P_area")
             continue
         end
 
@@ -504,7 +504,7 @@ function assemble!(problem::Problem{Mortar}, slave_element::Element{E}, time::Re
         # 6. add contribution to contact virtual work
         sdofs = get_gdofs(problem, slave_element)
         mdofs = get_gdofs(problem, master_element)
-        
+
         for i=1:field_dim
             lsdofs = sdofs[i:field_dim:end]
             lmdofs = mdofs[i:field_dim:end]
@@ -534,12 +534,12 @@ References
 
 """
 function assemble!(problem::Problem{Mortar}, slave_element::Element{E}, time::Real; first_slave_element=false) where E<:Union{Tri6}
-    
+
     props = problem.properties
     field_dim = get_unknown_field_dimension(problem)
     field_name = get_parent_field_name(problem)
     area = 0.0
-        
+
     Xs = slave_element("geometry", time)
 
     alp = props.alpha
@@ -554,7 +554,7 @@ function assemble!(problem::Problem{Mortar}, slave_element::Element{E}, time::Re
                 alp 0.0 alp 0.0 0.0 1.0-2*alp
             ]
     else
-        T = eye(6)
+        T = Matrix(1.0I, 6, 6)
     end
 
     #=
@@ -569,11 +569,11 @@ function assemble!(problem::Problem{Mortar}, slave_element::Element{E}, time::Re
     =#
 
     if props.dual_basis
-        # info("Creating dual basis for element $(slave_element.id)")            
+        # @info("Creating dual basis for element $(slave_element.id)")
         nsl = length(slave_element)
         De = zeros(nsl, nsl)
         Me = zeros(nsl, nsl)
-    
+
         # split slave element to linear sub-elements and loop
         for sub_slave_element in split_quadratic_element(slave_element, time)
 
@@ -587,7 +587,7 @@ function assemble!(problem::Problem{Mortar}, slave_element::Element{E}, time::Re
             N = vec(get_basis(sub_slave_element, xi, time))
             x0 = interpolate(N, X1)
             n0 = interpolate(N, n1)
-        
+
             # project slave nodes to auxiliary plane
             S = Vector[project_vertex_to_auxiliary_plane(X1[i], x0, n0) for i=1:nsl]
 
@@ -597,7 +597,7 @@ function assemble!(problem::Problem{Mortar}, slave_element::Element{E}, time::Re
             for master_element in master_elements
 
                 Xm = master_element("geometry", time)
-                
+
                 if norm(mean(Xs) - mean(Xm)) > problem.properties.distval
                     continue
                 end
@@ -632,28 +632,28 @@ function assemble!(problem::Problem{Mortar}, slave_element::Element{E}, time::Re
                             detJ = virtual_element(ip, time, Val{:detJ})
                             w = ip.weight*detJ
                             N1 = vec(slave_element(xi_s, time)*T)
-                            De += w*diagm(N1)
+                            De += w*Matrix(Diagonal(N1))
                             Me += w*N1*N1'
                         end
 
                     end # integration cells done
- 
-                end # sub aster elements done
+
+                end # sub master elements done
 
             end # master elements done
-     
+
         end # sub slave elements done
-        
+
         Ae = De*inv(Me)
-        # info("Dual basis construction finished.")
-        # info("Slave element geometry = $Xs")
-        # info("De = $De")
-        # info("Me = $Me")
-        # info("Dual basis coefficient matrix: $Ae")
+        # @info("Dual basis construction finished.")
+        # @info("Slave element geometry = $Xs")
+        # @info("De = $De")
+        # @info("Me = $Me")
+        # @info("Dual basis coefficient matrix: $Ae")
 
     else
         nsl = length(slave_element)
-        Ae = eye(nsl)
+        Ae = Matrix(1.0I, nsl, nsl)
     end
 
     # split slave element to linear sub-elements and loop
@@ -669,7 +669,7 @@ function assemble!(problem::Problem{Mortar}, slave_element::Element{E}, time::Re
         N = vec(get_basis(sub_slave_element, xi, time))
         x0 = interpolate(N, X1)
         n0 = interpolate(N, n1)
-    
+
         # project slave nodes to auxiliary plane
         S = Vector[project_vertex_to_auxiliary_plane(X1[i], x0, n0) for i=1:nsl]
 
@@ -679,7 +679,7 @@ function assemble!(problem::Problem{Mortar}, slave_element::Element{E}, time::Re
         for master_element in master_elements
 
             Xm = master_element("geometry", time)
- 
+
             if norm(mean(Xs) - mean(Xm)) > problem.properties.distval
                 continue
             end
@@ -702,7 +702,7 @@ function assemble!(problem::Problem{Mortar}, slave_element::Element{E}, time::Re
                 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)
-                    warn("Polygon P has zero area: $P_area")
+                    @warn("Polygon P has zero area: $P_area")
                     continue
                 end
 
@@ -752,7 +752,7 @@ function assemble!(problem::Problem{Mortar}, slave_element::Element{E}, time::Re
                 # 6. add contribution to contact virtual work
                 sdofs = get_gdofs(problem, slave_element)
                 mdofs = get_gdofs(problem, master_element)
-                
+
                 for i=1:field_dim
                     lsdofs = sdofs[i:field_dim:end]
                     lmdofs = mdofs[i:field_dim:end]
@@ -762,11 +762,11 @@ function assemble!(problem::Problem{Mortar}, slave_element::Element{E}, time::Re
                     add!(problem.assembly.C2, lsdofs, lmdofs, -Me)
                 end
                 add!(problem.assembly.g, sdofs, ge)
-   
+
             end # sub aster elements done
 
         end # master elements done
- 
+
     end # sub slave elements done
 
     return area
@@ -784,7 +784,7 @@ function assemble!(problem::Problem{Mortar}, time::Real, ::Type{Val{2}}, ::Type{
     #=
     if props.split_quadratic_slave_elements
         if !props.linear_surface_elements
-            warn("Mortar3D: split_quadratic_surfaces = true and linear_surface_elements = false maybe have unexpected behavior")
+            @warn("Mortar3D: split_quadratic_surfaces = true and linear_surface_elements = false maybe have unexpected behavior")
         end
         slave_elements = split_quadratic_elements(slave_elements, time)
     end
@@ -805,7 +805,7 @@ function assemble!(problem::Problem{Mortar}, time::Real, ::Type{Val{2}}, ::Type{
         first_slave_element = false
 
     end # slave elements done, contact virtual work ready
- 
+
     C1 = sparse(problem.assembly.C1)
     C2 = sparse(problem.assembly.C2)
 
@@ -843,13 +843,13 @@ function assemble!(problem::Problem{Mortar}, time::Real, ::Type{Val{2}}, ::Type{
         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)
+        d[get_nonzero_rows(T)] .= 0.0
+        T += sparse(Diagonal(d))
+        invT += sparse(Diagonal(d))
         #invT2 = sparse(inv(full(T)))
-        #info("invT == invT2? ", invT == invT2)
+        #@info("invT == invT2? ", invT == invT2)
         #maxabsdiff = maximum(abs(invT - invT2))
-        #info("max diff = $maxabsdiff")
+        #@info("max diff = $maxabsdiff")
         C1 = C1*invT
         C2 = C2*invT
     end
@@ -862,4 +862,3 @@ function assemble!(problem::Problem{Mortar}, time::Real, ::Type{Val{2}}, ::Type{
     problem.assembly.C2 = C2
 
 end
-
diff --git a/src/solvers.jl b/src/solvers.jl
index 5c77f3a..2b58d37 100644
--- a/src/solvers.jl
+++ b/src/solvers.jl
@@ -76,8 +76,8 @@ function get_field_assembly(solver::Solver)
     M = sparse(M, N, N)
     K = sparse(K, N, N)
     if nnz(K) == 0
-        warn("Field assembly seems to be empty. Check that elements are ",
-             "pushed to problem and formulation is correct.")
+        @warn("Field assembly seems to be empty. Check that elements are ",
+              "pushed to problem and formulation is correct.")
     end
     Kg = sparse(Kg, N, N)
     f = sparse(f, N, 1)
@@ -98,26 +98,26 @@ function check_for_overconstrained_dofs(solver::Solver)
         overconstrained_dofs = intersect(constrained_dofs, new_constraints)
         all_overconstrained_dofs = union(all_overconstrained_dofs, overconstrained_dofs)
         if length(overconstrained_dofs) != 0
-            warn("problem is overconstrained, finding overconstrained dofs... ")
+            @warn("problem is overconstrained, finding overconstrained dofs... ")
             overdetermined = true
             for dof in overconstrained_dofs
                 for problem_ in boundary_problems
                     new_constraints_ = Set(problem_.assembly.C2.I)
                     new_constraints_ = setdiff(new_constraints_, problem_.assembly.removed_dofs)
                     if dof in new_constraints_
-                        warn("overconstrained dof $dof defined in problem $(problem_.name)")
+                        @warn("overconstrained dof $dof defined in problem $(problem_.name)")
                     end
                 end
-                warn("To solve overconstrained situation, remove dofs from problems so that it exists only in one.")
-                warn("To do this, use push! to add dofs to remove to problem.assembly.removed_dofs, e.g.")
-                warn("`push!(bc.assembly.removed_dofs, $dof`)")
+                @warn("To solve overconstrained situation, remove dofs from problems so that it exists only in one.")
+                @warn("To do this, use push! to add dofs to remove to problem.assembly.removed_dofs, e.g.")
+                @warn("`push!(bc.assembly.removed_dofs, $dof`)")
             end
         end
         constrained_dofs = union(constrained_dofs, new_constraints)
     end
     if overdetermined
-        warn("List of all overconstrained dofs:")
-        warn(sort(collect(all_overconstrained_dofs)))
+        @warn("List of all overconstrained dofs:")
+        @warn(sort(collect(all_overconstrained_dofs)))
         error("problem is overconstrained, not continuing to solution.")
     end
     return true
@@ -149,9 +149,9 @@ function get_boundary_assembly(solver::Solver, N)
         f_ = sparse(assembly.f, N, 1)
         g_ = sparse(assembly.g, N, 1)
         for dof in assembly.removed_dofs
-            info("$(problem.name): removing dof $dof from assembly")
-            C1_[dof,:] = 0.0
-            C2_[dof,:] = 0.0
+            @info("$(problem.name): removing dof $dof from assembly")
+            C1_[dof,:] .= 0.0
+            C2_[dof,:] .= 0.0
         end
         SparseArrays.dropzeros!(C1_)
         SparseArrays.dropzeros!(C2_)
@@ -160,9 +160,9 @@ function get_boundary_assembly(solver::Solver, N)
         new_constraints = get_nonzero_rows(C2_)
         overconstrained_dofs = intersect(already_constrained, new_constraints)
         if length(overconstrained_dofs) != 0
-            warn("overconstrained dofs $overconstrained_dofs")
-            warn("already constrained = $already_constrained")
-            warn("new constraints = $new_constraints")
+            @warn("overconstrained dofs $overconstrained_dofs")
+            @warn("already constrained = $already_constrained")
+            @warn("new constraints = $new_constraints")
             overconstrained_dofs = sort(overconstrained_dofs)
             error("overconstrained dofs, not solving problem.")
         end
@@ -194,17 +194,17 @@ function solve!(solver::Solver, K, C1, C2, D, f, g, u, la, ::Type{Val{1}})
     I = setdiff(A, B)
 
     if length(B) == 0
-        warn("No rows in C2, forget to set Dirichlet boundary conditions to model?")
+        @warn("No rows in C2, forget to set Dirichlet boundary conditions to model?")
     else
-        u[B] = lufact(C2[B,B2]) \ full(g[B])
+        u[B] = lu(C2[B,B2]) \ Vector(g[B])
     end
 
     # solve interior domain using LDLt factorization
-    F = ldltfact(K[I,I])
-    u[I] = F \ (f[I] - K[I,B]*u[B])
+    F = ldlt(K[I,I])
+    u[I] = F \ Vector(f[I] - K[I,B]*u[B])
 
     # solve lagrange multipliers
-    la[B] = lufact(C1[B2,B]) \ full(f[B] - K[B,I]*u[I] - K[B,B]*u[B])
+    la[B] = lu(C1[B2,B]) \ Vector(f[B] - K[B,I]*u[I] - K[B,B]*u[B])
 
     return true
 end
@@ -249,14 +249,16 @@ function solve!(solver::Solver, K, C1, C2, D, f, g, u, la, ::Type{Val{3}})
     b = [f; g]
 
     ndofs = size(K, 2)
-    nz = ones(2*ndofs)
-    nz[get_nonzero_rows(A)] = 0.0
-    A += spdiagm(nz)
+    nonzero_rows = zeros(2*ndofs)
+    for j in rowvals(A)
+        nonzero_rows[j] = 1.0
+    end
+    A += sparse(Diagonal(1.0 .- nonzero_rows))
 
-    x = lufact(A) \ full(b)
+    x = lu(A) \ Vector(b[:])
 
-    u[:] = x[1:ndofs]
-    la[:] = x[ndofs+1:end]
+    u[:] .= x[1:ndofs]
+    la[:] .= x[ndofs+1:end]
 
     return true
 end
@@ -264,7 +266,7 @@ end
 """ Default linear system solver for solver. """
 function solve!(solver::Solver; empty_assemblies_before_solution=true, symmetric=true)
 
-    info("Solving problems ...")
+    @info("Solving linear system.")
     t0 = Base.time()
 
     # assemble field & boundary problems
@@ -288,7 +290,6 @@ function solve!(solver::Solver; empty_assemblies_before_solution=true, symmetric
         for problem in get_field_problems(solver)
             empty!(problem.assembly)
         end
-        gc()
     end
 
     #=
@@ -313,25 +314,25 @@ function solve!(solver::Solver; empty_assemblies_before_solution=true, symmetric
     u = zeros(ndofs)
     la = zeros(ndofs)
     is_solved = false
-    i = 0
+    local i
     for i in [1, 2, 3]
         is_solved = solve!(solver, K, C1, C2, D, f, g, u, la, Val{i})
         if is_solved
+            t1 = round(Base.time()-t0; digits=2)
+            norms = (norm(u), norm(la))
+            @info("Solved linear system in $t1 seconds using solver $i. " *
+                  "Solution norms (||u||, ||la||): $norms.")
             break
         end
     end
     if !is_solved
         error("Failed to solve linear system!")
     end
-    t1 = round(Base.time()-t0, 2)
-    norms = (norm(u), norm(la))
     #push!(solver.norms, norms)
-
     #solver.u = u
     #solver.la = la
 
-    info("Solved problems in $t1 seconds using solver $i.")
-    info("Solution norms = $norms.")
+    @info("")
 
     return u, la
 end
@@ -347,7 +348,7 @@ populated with global stiffness matrix, force vector, and, optionally,
 mass matrix.
 """
 function assemble!(solver::Solver, time::Float64; with_mass_matrix=false)
-    info("Assembling problems ...")
+    @info("Assembling problems ...")
 
     for problem in get_problems(solver)
         timeit("assemble $(problem.name)") do
@@ -373,7 +374,7 @@ function assemble!(solver::Solver, time::Float64; with_mass_matrix=false)
     end
     solver.ndofs = ndofs
     =#
-    info("Assembly done!")
+    @info("Assembly done!")
 end
 
 function get_unknown_fields(solver::Solver)
@@ -399,18 +400,18 @@ end
 """ Default initializer for solver. """
 function initialize!(solver::Solver)
     if solver.initialized
-        warn("initialize!(): solver already initialized")
+        @warn("initialize!(): solver already initialized")
         return
     end
-    info("Initializing solver ...")
+    @info("Initializing solver ...")
     problems = get_problems(solver)
     length(problems) != 0 || error("Empty solver, add problems to solver using push!")
     t0 = Base.time()
     field_problems = get_field_problems(solver)
-    length(field_problems) != 0 || warn("No field problem found from solver, add some..?")
+    length(field_problems) != 0 || @warn("No field problem found from solver, add some..?")
     field_name = get_unknown_field_name(solver)
     field_dim = get_unknown_field_dimension(solver)
-    info("initialize!(): looks we are solving $field_name, $field_dim dofs/node")
+    @info("initialize!(): looks we are solving $field_name, $field_dim dofs/node")
     nodes = Set{Int64}()
     for problem in problems
         initialize!(problem, solver.time)
@@ -420,9 +421,9 @@ function initialize!(solver::Solver)
         end
     end
     nnodes = length(nodes)
-    info("Total number of nodes in problems: $nnodes")
+    @info("Total number of nodes in problems: $nnodes")
     maxdof = maximum(nodes)*field_dim
-    info("# of max dof (=size of solution vector) is $maxdof")
+    @info("# of max dof (=size of solution vector) is $maxdof")
     solver.u = zeros(maxdof)
     solver.la = zeros(maxdof)
     # TODO: this could be used to initialize elements too...
@@ -432,8 +433,8 @@ function initialize!(solver::Solver)
         problem.assembly.la = zeros(maxdof)
         # initialize(problem, ....)
     end
-    t1 = round(Base.time()-t0, 2)
-    info("Initialized solver in $t1 seconds.")
+    t1 = round(Base.time()-t0; digits=2)
+    @info("Initialized solver in $t1 seconds.")
     solver.initialized = true
 end
 
@@ -459,7 +460,7 @@ function (solver::Solver)(field_name::String, time::Float64)
             continue
         end
         if length(field) == 0
-            warn("no field $field_name found for problem $(problem.name)")
+            @warn("no field $field_name found for problem $(problem.name)")
             continue
         end
         push!(fields, field)
@@ -470,14 +471,7 @@ function (solver::Solver)(field_name::String, time::Float64)
     return merge(fields...)
 end
 
-""" Default update for solver. """
 function update!(solver::Solver{S}, u, la, time) where S
-    #u = solver.u
-    #la = solver.la
-
-    info("Updating problems ...")
-    t0 = Base.time()
-
     for problem in get_problems(solver)
         assembly = get_assembly(problem)
         elements = get_elements(problem)
@@ -486,9 +480,6 @@ function update!(solver::Solver{S}, u, la, time) where S
         # .. and then from assembly (u,la) to elements
         update!(problem, assembly, elements, time)
     end
-
-    t1 = round(Base.time()-t0, 2)
-    info("Updated problems in $t1 seconds.")
 end
 
 """ Default postprocess for solver. Loop all problems and run postprocess
@@ -496,11 +487,13 @@ functions to calculate secondary fields, i.e. contact pressure, stress,
 heat flux, reaction force etc. quantities.
 """
 function postprocess!(solver::Solver, time)
-    info("Running postprocess scripts for solver...")
-    for problem in get_problems(solver)
+    problems = get_problems(solver)
+    nproblems = length(problems)
+    @info("Postprocessing $nproblems problems.")
+    for problem in problems
         for field_name in problem.postprocess_fields
             field = Val{Symbol(field_name)}
-            info("Running postprocess for problem $(problem.name), field $field_name")
+            @info("Running postprocess for problem $(problem.name), field $field_name")
             postprocess!(problem, time, field)
         end
     end
@@ -517,9 +510,10 @@ to Xdmf file. By default write the main unknown field (displacement, temperature
 function write_results!(solver, time)
     results_writers = get_results_writers(solver)
     if length(results_writers) == 0
-        info("Xdmf is not attached to solver, not writing output to a file.")
-        info("To write results to Xdmf file, attach Xdmf to Solver, i.e.")
-        info("add_results_writer!(solver, Xdmf(\"results\"))")
+        @info("No result writers are attached to analysis, not writing output.")
+        @info("To write results to Xdmf file, attach Xdmf to analysis, i.e.")
+        @info("xdmf_output = Xdmf(\"simulation_results\")")
+        @info("add_results_writer!(analysis, xdmf_output)")
         return
     end
     # FIXME: result writer can be anything, not only Xdmf
@@ -585,10 +579,10 @@ function FEMBase.run!(solver::Solver{Nonlinear})
 
     # 2. start non-linear iterations
     for properties.iteration=1:properties.max_iterations
-        info(repeat("-", 80))
-        info("Starting nonlinear iteration #$(properties.iteration)")
-        info("Increment time t=$(round(time, 3))")
-        info(repeat("-", 80))
+        @info(repeat("-", 80))
+        @info("Starting nonlinear iteration #$(properties.iteration)")
+        @info("Increment time t=$(round(time; digits=3))")
+        @info(repeat("-", 80))
 
         # 2.1 update assemblies
         for problem in problems
@@ -604,7 +598,7 @@ function FEMBase.run!(solver::Solver{Nonlinear})
 
         # 2.4 check convergence
         if properties.iteration >= properties.min_iterations && has_converged(solver)
-            info("Converged in $(properties.iteration) iterations.")
+            @info("Converged in $(properties.iteration) iterations.")
             # 2.4.1 run any postprocessing of problems
             postprocess!(solver, time)
             # 2.4.2 update Xdmf output
@@ -638,17 +632,22 @@ function Linear()
     return Linear(0.0)
 end
 
-function FEMBase.run!(solver::Analysis{Linear})
-    time = solver.properties.time
-    problems = get_problems(solver)
-    N = 0
+function FEMBase.run!(analysis::Analysis{Linear})
+    time = analysis.properties.time
+    @info("Running linear quasistatic analysis `$(analysis.name)` at time $time.")
+    problems = get_problems(analysis)
+    nproblems = length(problems)
+    @info("Assembling $nproblems problems.")
     @timeit "assemble problems" for problem in problems
         isempty(problem.assembly) || continue
         initialize!(problem, time)
         assemble!(problem, time)
     end
-    @timeit "solve linear system" u, la = solve!(solver)
-    @timeit "update problems" update!(solver, u, la, time)
+    @timeit "solve linear system" u, la = solve!(analysis)
+    @timeit "update problems" update!(analysis, u, la, time)
+    postprocess!(analysis, time)
+    write_results!(analysis, time)
+    @info("Quasistatic linear analysis ready.")
 end
 
 # Convenience functions
@@ -676,13 +675,13 @@ end
 # will be deprecated
 
 function (solver::Solver)(time::Float64=0.0)
-    warn("analysis(time) is deprecated. Instead, use run!(analysis)")
+    @warn("analysis(time) is deprecated. Instead, use run!(analysis)")
     solver.properties.time = time
     run!(solver)
 end
 
 function solve!(solver::Solver, time::Float64)
-    warn("solve!(analysis, time) is deprecated. Instead, use run!(analysis)")
+    @warn("solve!(analysis, time) is deprecated. Instead, use run!(analysis)")
     solver.properties.time = time
     run!(solver)
 end
diff --git a/src/solvers_modal.jl b/src/solvers_modal.jl
index b465277..4ac9038 100644
--- a/src/solvers_modal.jl
+++ b/src/solvers_modal.jl
@@ -1,16 +1,8 @@
 # This file is a part of JuliaFEM.
 # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
 
-""" Modal solver to solve generalized eigenvalue problems Ku = Muλ
+using SparseArrays, Arpack
 
-Examples
---------
-
-julia> problems = get_problems()
-julia> solver = Solver(Modal)
-julia> push!(solver, problems...)
-julia> solver()
-"""
 mutable struct Modal <: AbstractSolver
     time :: Float64
     geometric_stiffness :: Bool
@@ -28,7 +20,7 @@ mutable struct Modal <: AbstractSolver
 end
 
 function Modal(nev=10, which=:SM)
-    solver = Modal(0.0, false, [], Matrix{Float64}(0,0), nev, which,
+    solver = Modal(0.0, false, [], Matrix{Float64}(undef,0,0), nev, which,
                    false, [], true, true, false, false, 0.0)
 end
 
@@ -53,8 +45,8 @@ function calc_projection(problem::T) where
     M = -C2[s,m]
 
     if !isdiag(D)
-        warn("Mortar matrix D is not diagonal. This might take a long time.")
-        P = ldltfact(1/2*(D + D')) \ M
+        @warn("Mortar matrix D is not diagonal. This might take a long time.")
+        P = ldlt(1/2*(D + D')) \ M
     else
         P = D \ M
     end
@@ -68,10 +60,12 @@ function FEMBase.eliminate_boundary_conditions!(problem::P, K, M, f) where {P}
     C2 = sparse(problem.assembly.C2)
     C1 == C2 || error("Cannot eliminate boundary condition $P: C1 != C2.")
     isdiag(C1) || error("Cannot eliminate boundary condition $P: C is not diagonal")
-    info("Eliminating boundary condition $(problem.name) from global system.")
+    @info("Eliminating boundary condition $(problem.name) from global system.")
     fixed_dofs = get_nonzero_rows(C1)
-    K[fixed_dofs,:] = K[:,fixed_dofs] = 0.0
-    M[fixed_dofs,:] = M[:,fixed_dofs] = 0.0
+    K[fixed_dofs,:] .= 0.0
+    K[:,fixed_dofs] .= 0.0
+    M[fixed_dofs,:] .= 0.0
+    M[:,fixed_dofs] .= 0.0
     dropzeros!(K)
     dropzeros!(M)
     return nothing
@@ -84,15 +78,15 @@ Eliminate Mortar boundary condition from matrices K, M and force vector f.
 """
 function FEMBase.eliminate_boundary_conditions!(problem::T, K, M, f) where
                             {T <: Union{Problem{Mortar}, Problem{Mortar2D}}}
-    info("Eliminating mesh tie constraint $(problem.name) using static condensation")
+    @info("Eliminating mesh tie constraint $(problem.name) using static condensation")
     s, m, P = calc_projection(problem)
     ndim = size(K, 1)
     Id = ones(ndim)
-    Id[s] = 0.0
-    Q = spdiagm(Id)
+    Id[s] .= 0.0
+    Q = sparse(Diagonal(Id))
     Q[s,m] += P
-    K[:,:] = Q'*K*Q
-    M[:,:] = Q'*M*Q
+    K[:,:] .= Q'*K*Q
+    M[:,:] .= Q'*M*Q
     return nothing
 end
 
@@ -100,10 +94,8 @@ function FEMBase.run!(solver::Solver{Modal})
     time = solver.properties.time
     problems = get_problems(solver)
     properties = solver.properties
-    info(repeat("-", 80))
-    info("Starting natural frequency solver")
-    info("Increment time t=$(round(time, 3))")
-    info(repeat("-", 80))
+
+    @info("Starting natural frequency solver at time $time")
 
     @timeit "assemble matrices" begin
         assemble!(solver, time; with_mass_matrix=true)
@@ -116,17 +108,14 @@ function FEMBase.run!(solver::Solver{Modal})
     dim = size(K, 1)
     ndofs = size(K, 1)
 
-    K_red = K
-    M_red = M
-
     for P in properties.P
-        info("Using P to make transformation K_red = P'*K*P and M_red = P'*M*P")
-        K_red[:,:] = P'*K_red*P
-        M_red[:,:] = P'*M_red*P
+        @info("Using P to make transformation K_red = P'*K*P and M_red = P'*M*P")
+        K[:,:] .= P'*K*P
+        M[:,:] .= P'*M*P
     end
 
     for problem in get_problems(solver)
-        eliminate_boundary_conditions!(problem, K_red, M_red, f)
+        eliminate_boundary_conditions!(problem, K, M, f)
     end
 
     # free up some memory before solution
@@ -134,42 +123,39 @@ function FEMBase.run!(solver::Solver{Modal})
         for problem in get_field_problems(solver)
             empty!(problem.assembly)
         end
-        gc()
     end
 
-    SparseArrays.droptol!(K_red, 1.0e-9)
-    SparseArrays.droptol!(M_red, 1.0e-9)
-    nz = get_nonzero_rows(K_red)
-    K_red = K_red[nz,nz]
-    M_red = M_red[nz,nz]
+    SparseArrays.droptol!(K, 1.0e-9)
+    SparseArrays.droptol!(M, 1.0e-9)
+    nz = get_nonzero_rows(K)
+    K = K[nz,nz]
+    M = M[nz,nz]
 
     sigma = 0.0
     if properties.sigma != 0.0
-        info("Adding diagonal term $(properties.sigma) to stiffness matrix")
+        @info("Adding diagonal term $(properties.sigma) to stiffness matrix")
         sigma = properties.sigma
     end
 
     props = solver.properties
 
-    info("Calculate $(props.nev) eigenvalues...")
-
-    tic()
+    @debug("Calculating $(props.nev) eigenvalues...")
 
     if properties.symmetric
-        K_red = 1/2*(K_red + transpose(K_red))
-        M_red = 1/2*(M_red + transpose(M_red))
+        K = Symmetric(K)
+        M = Symmetric(M)
     end
 
     if properties.info_matrices
-        info("is K symmetric? ", issymmetric(K_red))
-        info("is M symmetric? ", issymmetric(M_red))
-        info("is K positive definite? ", isposdef(K_red))
-        info("is M positive definite? ", isposdef(M_red))
+        @info("is K symmetric? ", issymmetric(K))
+        @info("is M symmetric? ", issymmetric(M))
+        @info("is K positive definite? ", isposdef(K))
+        @info("is M positive definite? ", isposdef(M))
     end
 
     if properties.dense
-        K_red = full(K_red)
-        M_red = full(M_red)
+        K = Matrix(K)
+        M = Matrix(M)
     end
 
     om2 = nothing
@@ -178,49 +164,39 @@ function FEMBase.run!(solver::Solver{Modal})
 
     try
         @timeit "solve eigenvalue problem using `eigs`" begin
-            om2, X = eigs(K_red + sigma*I, M_red; nev=props.nev, which=props.which)
+            om2, X = eigs(K + sigma*I, M; nev=props.nev, which=props.which)
         end
         passed = true
     catch
-        info("Failed to calculate eigenvalues for problem.")
-        b1 = issymmetric(K_red)
-        b2 = issymmetric(M_red)
-        b3 = isposdef(K_red)
-        b4 = isposdef(M_red)
-        info("Is K symmetric? $b1")
-        info("Is M symmetric? $b2")
-        info("Is K positive definite? $b3")
-        info("Is M positive definite? $b4")
-        if properties.sigma != 0.0
-            if !b3
-                info("Stiffness matrix is not positive definite and Cholesky ",
-                     "factorization is failing. Model is not supported enough ",
-                     "with boundary conditions. To work around this problem, ",
-                     "use `problem.properties.sigma = <some small value>` ",
-                     "To add artificial stiffness to model. (Or add boundary ",
-                     "conditions.)")
-            end
+        @info("Failed to calculate eigenvalues for problem.",
+              issymmetric(K), issymmetric(M), isposdef(K), isposdef(M))
+        if !isapprox(properties.sigma, 0.0)
+            @info("Stiffness matrix is not positive definite and Cholesky " *
+                  "factorization is failing. Model is not supported enough " *
+                  "with boundary conditions. To work around this problem, " *
+                  "use `problem.properties.sigma = <some small value>` " *
+                  "To add artificial stiffness to model. (Or add boundary " *
+                  "conditions.)")
             rethrow()
         end
     end
 
     if !passed
         sigma = props.sigma = 1.0e-9
-        info("Calculation of eigenvalues failed. Stiffness matrix is not ",
-             "positive definite and Cholesky factorization is failing. Trying ",
-             "again by adjusting problem.properties.sigma to $sigma.")
+        @info("Calculation of eigenvalues failed. Stiffness matrix is not " *
+              "positive definite and Cholesky factorization is failing. Trying " *
+              "again by adjusting problem.properties.sigma to $sigma.")
         try
-            om2, X = eigs(K_red + sigma*I, M_red; nev=props.nev, which=props.which)
+            om2, X = eigs(K + sigma*I, M; nev=props.nev, which=props.which)
             passed = true
         catch
-            info("Failed to calculate eigenvalues with sigma value $sigma. ",
-                 "Manually set sigma to something larger and try again.")
+            @info("Failed to calculate eigenvalues with sigma value $sigma. " *
+                  "Manually set sigma to something larger and try again.")
             rethrow()
         end
     end
 
-    t1 = round(toq(), 2)
-    info("Eigenvalues computed in $t1 seconds. Squared eigenvalues: $om2")
+    @info("Squared eigenvalues: $om2.")
 
     props.eigvals = om2
     neigvals = length(om2)
@@ -239,7 +215,7 @@ function FEMBase.run!(solver::Solver{Modal})
 
     @timeit "save results to Xdmf" update_xdmf!(solver)
 
-    return true
+    return nothing
 
 end
 
@@ -247,13 +223,13 @@ function update_xdmf!(solver::Solver{Modal})
 
     results_writers = get_results_writers(solver)
     if length(results_writers) == 0
-        info("Xdmf is not attached to solver, not writing output to a file.")
-        info("To write results to Xdmf file, attach Xdmf to Solver, i.e.")
-        info("add_results_writer!(solver, Xdmf(\"results\"))")
+        @info("Xdmf is not attached to solver, not writing output to a file.")
+        @info("To write results to Xdmf file, attach Xdmf to Solver, i.e.")
+        @info("add_results_writer!(solver, Xdmf(\"results\"))")
         return
     end
     if maximum(abs.(imag(solver.properties.eigvals))) > 1.0e-9
-        info("Writing imaginary eigenvalues for Xdmf not supported.")
+        @info("Writing imaginary eigenvalues for Xdmf not supported.")
         return
     end
 
@@ -265,9 +241,9 @@ function update_xdmf!(solver::Solver{Modal})
     nnodes = length(X_)
     ndofs = round(Int, size(solver.properties.eigvecs, 1)/nnodes)
     ndim = length(X_[first(node_ids)])
-    info("Number of nodes: $nnodes. ",
-         "Number of dofs/node: $ndofs. ",
-         "Dimension of geometry: $ndim.")
+    @info("Number of nodes: $nnodes. ",
+          "Number of dofs/node: $ndofs. ",
+          "Dimension of geometry: $ndim.")
     @timeit "create ncoords array" begin
         X = zeros(ndim, nnodes)
         for j in node_ids
@@ -318,12 +294,12 @@ function update_xdmf!(solver::Solver{Modal})
 
     @timeit "save modes" for (j, eigval) in enumerate(real(solver.properties.eigvals))
         if eigval < 0.0
-            warn("negative real eigenvalue found, om2=$eigval, setting to zero.")
+            @warn("negative real eigenvalue found, om2=$eigval, setting to zero.")
             eigval = 0.0
         end
         freq = sqrt(eigval)/(2.0*pi)
         path = "/Results/Natural Frequency Analysis/$unknown_field_name/Mode $j"
-        info("Creating frequency frame f=$(round(freq, 3)), path=$path")
+        @info("Creating frequency frame f=$(round(freq; digits=3)), path=$path")
 
         frame = new_element("Grid")
         time = new_child(frame, "Time")
@@ -393,7 +369,7 @@ function update_xdmf!(solver::Solver{Modal})
 end
 
 function solve!(solver::Solver{Modal}, time::Float64)
-    info("solve!(analysis, time) is deprecated. Use run!(analysis) instead.")
+    @info("solve!(analysis, time) is deprecated. Use run!(analysis) instead.")
     solver.properties.time = time
     run!(solver)
 end