Fix deprecation warnings

* Add docstrings * Refactor code * Module level docstring giving an example
2026-09-10 04:50:49 +00:00 · 2018-09-06 12:02:56 +03:00
parent 1357e7e46d
commit ca911e5e87
15 changed files with 539 additions and 432 deletions
@@ -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")

@@ -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
@@ -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
@@ -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
@@ -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})
@@ -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


@@ -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
@@ -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

@@ -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)
@@ -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

@@ -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]
@@ -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
@@ -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
-
@@ -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
@@ -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