new style dict field, xdmf improvements

src/JuliaFEM.jl

@@ -1,6 +1,12 @@
 # This file is a part of JuliaFEM.
 # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
 
+using Logging
+
+if haskey(ENV, "JULIAFEM_LOGLEVEL")
+    ENV["JULIAFEM_LOGLEVEL"] == "DEBUG" && Logging.configure(level=DEBUG)
+end
+
 """
 This is JuliaFEM -- Finite Element Package
 """

src/elements.jl

@@ -11,10 +11,12 @@ type Element{E<:AbstractElement}
     properties :: E
 end
 
-function Element{E<:AbstractElement}(::Type{E}, connectivity=[], integration_points=[], id=-1, fields=Dict(), properties...)
-    variant = E(properties...)
-    element = Element{E}(id, connectivity, integration_points, fields, variant)
-    return element
+function Element{E<:AbstractElement}(::Type{E}, id::Int64, connectivity::Vector{Int64})
+    return Element{E}(id, connectivity, [], Dict(), E())
+end
+
+function Element{E<:AbstractElement}(::Type{E}, connectivity::Vector{Int64})
+    return Element{E}(-1, connectivity, [], Dict(), E())
 end
 
 function getindex(element::Element, field_name::AbstractString)
@@ -59,9 +61,15 @@ function call(element::Element, ip, time::Float64=0.0)
 end
 
 function call(element::Element, ip, time::Float64, ::Type{Val{:Jacobian}})
-    X = element["geometry"](time)
+    X = element("geometry", time)
     dN = get_dbasis(element, ip, time)
-    J = sum([kron(dN[:,i], X[i]') for i=1:length(X)])
+    nbasis = length(element)
+    if isa(X.data, Vector)
+        J = sum([kron(dN[:,i], X[i]') for i=1:nbasis])
+    else
+        c = get_connectivity(element)
+        J = sum([kron(dN[:,i], X[c[i]]') for i=1:nbasis])
+    end
     return J
 end
 
@@ -122,11 +130,16 @@ function call(element::Element, field::Field, ip, time::Float64)
     field_ = field(time)
     basis = element(ip, time)
     n = length(element)
-    m = length(field_)
-    if n != m
-        error("Error when trying to interpolate field $field at coords $ip and time $time: element length is $n and field length is $m, f = Nᵢfᵢ makes no sense!")
+    if isa(field_.data, Vector)
+        m = length(field_)
+        if n != m
+            error("Error when trying to interpolate field $field at coords $ip and time $time: element length is $n and field length is $m, f = Nᵢfᵢ makes no sense!")
+        end
+        return sum([field_[i]*basis[i] for i=1:n])
+    else
+        c = get_connectivity(element)
+        return sum([field_[c[i]]*basis[i] for i=1:n])
     end
-    return sum([field_[i]*basis[i] for i=1:n])
 end
 
 function size(element::Element, dim)
@@ -145,13 +158,19 @@ As a result element now have time invariant (variable) vector field "geometry" w
 
 """
 function update!(element::Element, field_name, data::Dict)
-    element[field_name] = [data[i] for i in get_connectivity(element)]
+    element[field_name] = Field(data)
+    #element[field_name] = [data[i] for i in get_connectivity(element)]
 end
 
 function update!{K,V}(element::Element, field_name, data::Pair{Float64, Dict{K, V}})
-    time, field_data = data
-    element_data = V[field_data[i] for i in get_connectivity(element)]
-    update!(element, field_name, time => element_data)
+    #time, field_data = data
+    #element_data = V[field_data[i] for i in get_connectivity(element)]
+    #update!(element, field_name, time => element_data)
+    if haskey(element, field_name)
+        update!(element[field_name], data)
+    else
+        element[field_name] = Field(data)
+    end
 end
 
 function update!(element::Element, field_name::AbstractString, datas::Union{Real, Vector, Pair{Float64, Union{Float64, Real, Vector{Any}}}}...)
@@ -321,4 +340,3 @@ function inside{E}(element::Element{E}, X, time)
     xi = get_local_coordinates(element, X, time)
     return inside(E, xi)
 end
-

src/fields.jl

@@ -97,11 +97,24 @@ end
 function Field{T}(data::Pair{Float64, T}...)
     return DCTV([Increment{T}(d[1], d[2]) for d in data])
 end
-
+#=
 function Field{T}(data::Pair{Float64, Vector{T}}...)
     return DVTV([Increment{Vector{T}}(d[1], d[2]) for d in data])
 end
 
+function Field{T}(data::Pair{Float64, Dict{Int64, T}}...)
+    return DVTV([Increment{Dict{Int64, T}}(d[1], d[2]) for d in data])
+end
+=#
+
+function Field{T<:Union{Vector, Dict}}(data::Pair{Float64, T}...)
+    return DVTV([Increment{T}(d[1], d[2]) for d in data])
+end
+
+function Field(data::Dict)
+    return DVTI(data)
+end
+
 function convert{T}(::Type{DCTV}, data::Pair{Real, Vector{T}}...)
     return DCTV([Increment{Vector{T}}(d[1], d[2]) for d in data])
 end
@@ -222,11 +235,11 @@ function Base.(:*)(N::Matrix, f::DCTI)
 end
 
 
-#  
+#
 #   Multiply DVTI field with another vector T. Vector length
 #   must match to the field length and this can be used mainly
 #   for interpolation purposes, i.e., u = ∑ Nᵢuᵢ
-#  
+#
 function Base.(:*)(T::Vector, f::DVTI)
     @assert length(T) == length(f)
     return sum([T[i]*f[i] for i=1:length(f)])
@@ -430,3 +443,7 @@ end
 function keys(field::DVTV)
     return Float64[increment.time for increment in field]
 end
+
+function setindex!(field::Field, val, idx::Int64)
+    field.data[idx] = val
+end

src/io.jl

@@ -22,6 +22,42 @@ function haskey(x::XMLElement, key::AbstractString)
     return has_child(x, key) || has_attribute(x, key)
 end
 
+function has_child(x::XMLElement, child_name::AbstractString)
+    return get_child(x, child_name) != nothing
+end
+
+function get_attribute(x::XMLElement, attr_name::AbstractString)
+    attr = attribute(x, attr_name)
+    numeric = tryparse(Int64, attr)
+    isnull(numeric) && (numeric = tryparse(Float64, attr))
+    isnull(numeric) && return attr
+    return get(numeric)
+end
+
+function new_child(xparent::XMLElement, name::AbstractString, attrs::Dict)
+    x = new_child(xparent, name)
+    for (k, v) in attrs
+        x[k] = v
+    end
+    return x
+end
+
+function new_child(xparent::XMLElement, name::AbstractString, attrs::Pair...)
+    x = new_child(xparent, name)
+    for (k, v) in attrs
+        x[k] = v
+    end
+    return x
+end
+
+""" Basic traverse support, so that it's possible to find data from xml using
+path syntax e.g. /foo/bar[2]/baz[@Name=Frame 1]/DataItem. If several elements
+with same name exists in tree, pick first by default and next ones can be picked
+using [] syntax or [@attr=value] syntax, see [1] for details. For last item use
+[end].
+
+[1] http://www.xdmf.org/index.php/XDMF_Model_and_Format
+"""
 function get_child(x::XMLElement, child_name::AbstractString)
     '/' in child_name && return nothing
     m = match(r"(\w+)\[(.+)\]", child_name)
@@ -52,18 +88,6 @@ function get_child(x::XMLElement, child_name::AbstractString)
     throw("Unable to parse: $(m[2])")
 end
 
-function has_child(x::XMLElement, child_name::AbstractString)
-    return get_child(x, child_name) != nothing
-end
-
-function get_attribute(x::XMLElement, attr_name::AbstractString)
-    attr = attribute(x, attr_name)
-    numeric = tryparse(Int64, attr)
-    isnull(numeric) && (numeric = tryparse(Float64, attr))
-    isnull(numeric) && return attr
-    return get(numeric)
-end
-
 function getindex(x::XMLElement, attr_name::AbstractString)
     attr_name = strip(attr_name, '/')
     child = get_child(x, attr_name)
@@ -82,22 +106,6 @@ function getindex(x::XMLElement, attr_name::AbstractString)
     end
 end
 
-function new_child(xparent::XMLElement, name::AbstractString, attrs::Dict)
-    x = new_child(xparent, name)
-    for (k, v) in attrs
-        x[k] = v
-    end
-    return x
-end
-
-function new_child(xparent::XMLElement, name::AbstractString, attrs::Pair...)
-    x = new_child(xparent, name)
-    for (k, v) in attrs
-        x[k] = v
-    end
-    return x
-end
-
 type Xdmf
     name :: AbstractString
     xml :: XMLElement

src/postprocess_utils.jl

@@ -199,28 +199,6 @@ function copy_field!(src_problem::Problem, dst_problem::Problem, field_name, tim
     copy_field!(src_problem.elements, dst_problem.elements, field_name, time)
 end
 
-""" Return field calculated to nodal points for elements in problem p. """
-function call(problem::Problem, field_name::AbstractString, time::Float64=0.0)
-    f = nothing
-    for element in get_elements(problem)
-        haskey(element, field_name) || continue
-        for (c, v) in zip(get_connectivity(element), element(field_name, time))
-            if f == nothing
-                f = Dict(c => v)
-                continue
-            end
-            if haskey(f, c)
-                if !isapprox(f[c], v)
-                    info("several values for single node when returning field $field_name")
-                    info("already have: $(f[c]), and trying to set $v")
-                end
-            else
-                f[c] = v
-            end
-        end
-    end
-    return f
-end
 
 function to_dataframe(u::Dict, abbreviation::Symbol)
     length(u) != 0 || return DataFrame()
@@ -419,20 +397,3 @@ function calculate_second_moment_of_mass(problem::Problem, X=[0.0, 0.0, 0.0], ti
     end
     return I
 end
-
-function getindex(problem::Problem, field_name::AbstractString)
-    info("fetching result $field_name")
-    timeframes = []
-    for frame in first(problem.elements)[field_name].data
-        push!(timeframes, frame.time)
-    end
-    info("time frames: $timeframes")
-    conn = get_connectivity(problem)
-    increments = Increment[]
-    for time in timeframes
-        p = problem(field_name, time)
-        data = [p[id] for id in conn]
-        push!(increments, Increment(time, data))
-    end
-    return DVTV(increments)
-end

src/problems.jl

@@ -8,8 +8,8 @@ abstract MixedProblem <: AbstractProblem
 
 """
 General linearized problem to solve
-    (K₁+K₂)*Δu +   C1.T*λ = f₁+f₂
-         C2*Δu +      D*λ = g
+    (K₁+K₂)Δu  +   C1*Δλ = f₁+f₂
+         C2Δu  +    D*Δλ = g
 """
 type Assembly
 
@@ -19,7 +19,7 @@ type Assembly
     K :: SparseMatrixCOO   # stiffness matrix
     Kg :: SparseMatrixCOO  # geometric stiffness matrix
     f :: SparseMatrixCOO   # force vector
-    fg :: SparseMatrixCOO  # 
+    fg :: SparseMatrixCOO  #
 
     # for boundary assembly
     C1 :: SparseMatrixCOO
@@ -90,6 +90,7 @@ type Problem{P<:AbstractProblem}
     elements :: Vector{Element}
     dofmap :: Dict{Element, Vector{Int64}} # connects element local dofs to global dofs
     assembly :: Assembly
+    fields :: Dict{AbstractString, Field}
     properties :: P
 end
 
@@ -104,10 +105,10 @@ julia> prob2 = Problem(Elasticity, 3)
 
 """
 function Problem{P<:FieldProblem}(::Type{P}, name::AbstractString, dimension::Int64)
-    return Problem{P}(name, dimension, "none", [], Dict(), Assembly(), P())
+    return Problem{P}(name, dimension, "none", [], Dict(), Assembly(), Dict(), P())
 end
 function Problem{P<:FieldProblem}(::Type{P}, dimension::Int64)
-    return Problem{P}("$P problem", dimension, "none", [], Dict(), Assembly(), P())
+    return Problem{P}("$P problem", dimension, "none", [], Dict(), Assembly(), Dict(), P())
 end
 
 """ Construct a new boundary problem.
@@ -117,16 +118,16 @@ Examples
 Create Dirichlet boundary problem for vector-valued (dim=3) elasticity problem.
 
 julia> bc1 = Problem(Dirichlet, "support", 3, "displacement")
-
+solver.
 """
 function Problem{P<:BoundaryProblem}(::Type{P}, name, dimension, parent_field_name)
-    return Problem{P}(name, dimension, parent_field_name, [], Dict(), Assembly(), P())
+    return Problem{P}(name, dimension, parent_field_name, [], Dict(), Assembly(), Dict(), P())
 end
 function Problem{P<:BoundaryProblem}(::Type{P}, main_problem::Problem)
     name = "$P problem"
     dimension = get_unknown_field_dimension(main_problem)
     parent_field_name = get_unknown_field_name(main_problem)
-    return Problem{P}(name, dimension, parent_field_name, [], Dict(), Assembly(), P())
+    return Problem{P}(name, dimension, parent_field_name, [], Dict(), Assembly(), Dict(), P())
 end
 
 function get_formulation_type{P<:FieldProblem}(problem::Problem{P})
@@ -198,6 +199,7 @@ function update!(problem::Problem, assembly::Assembly, u::Vector, la::Vector; ve
     # incremental formulation we solve KΔu = f and u = u + Δu
     assembly.u_prev = copy(assembly.u)
     assembly.la_prev = copy(assembly.la)
+
     if get_formulation_type(problem) == :total
         verbose && info("$(problem.name): total formulation, replacing solution vector with new values")
         assembly.u = u
@@ -225,7 +227,7 @@ end
 
 Notes
 -----
-If length of solution vector != number of nodes, i.e. field dimension is 
+If length of solution vector != number of nodes, i.e. field dimension is
 something other than 1, reshape vectors so it's length matches to the
 number of nodes so that one can easily get nodal results.
 """
@@ -282,9 +284,50 @@ function length(problem::Problem)
 end
 
 function update!(problem::Problem, field_name::AbstractString, data)
+    if haskey(problem.fields, field_name)
+        update!(problem.fields[field_name], field_name::AbstractString, data)
+    else
+        problem.fields[field_name] = Field(data)
+    end
     update!(problem.elements, field_name::AbstractString, data)
 end
 
+function haskey(problem::Problem, field_name::AbstractString)
+    return haskey(problem.fields, field_name)
+end
+
+function getindex(problem::Problem, field_name::AbstractString)
+    return problem.fields[field_name]
+end
+
+""" Return field calculated to nodal points for elements in problem p. """
+function call(problem::Problem, field_name::AbstractString, time::Float64=0.0)
+    if haskey(problem, field_name)
+        return problem[field_name](time)
+    end
+    f = nothing
+    for element in get_elements(problem)
+        haskey(element, field_name) || continue
+        for (c, v) in zip(get_connectivity(element), element(field_name, time))
+            if f == nothing
+                f = Dict(c => v)
+                continue
+            end
+            if haskey(f, c)
+                if !isapprox(f[c], v)
+                    info("several values for single node when returning field $field_name")
+                    info("already have: $(f[c]), and trying to set $v")
+                end
+            else
+                f[c] = v
+            end
+        end
+    end
+    f == nothing && return f
+    update!(problem, field_name, time => f)
+    return f
+end
+
 """ Return the dimension of the unknown field of this problem. """
 function get_unknown_field_dimension(problem::Problem)
     return problem.dimension
@@ -381,4 +424,3 @@ function find_nodes_by_dofs(dim, dofs)
     end
     return nodes
 end
-

src/problems_elasticity.jl

@@ -4,7 +4,7 @@
 """ Elasticity equations.
 
 Field equation is:
-    
+
     m∂²u/∂t² = ∇⋅σ - b
 
 Weak form is: find u∈U such that ∀v in V
@@ -151,7 +151,7 @@ function assemble{El<:Elasticity2DVolumeElements}(problem::Problem{Elasticity},
         # cauchy_stress = F'*stress*F/det(F)
         # cauchy_stress = [cauchy_stress[1,1]; cauchy_stress[2,2]; cauchy_stress[1,2]]
         # update!(ip, "cauchy stress", time => cauchy_stress)
-        
+
         # material stiffness end
 
         if props.geometric_stiffness
@@ -170,13 +170,13 @@ function assemble{El<:Elasticity2DVolumeElements}(problem::Problem{Elasticity},
             S2[2,2] = stress_vec[2]
             S2[1,2] = S2[2,1] = stress_vec[3]
             S2[3:4,3:4] = S2[1:2,1:2]
-            
+
             Kg += w*BNL'*S2*BNL # geometric stiffness
 
         end
 
         # rhs, internal and external load
-        
+
         f -= w*BL'*stress_vec
 
         if haskey(element, "displacement load")
@@ -747,4 +747,3 @@ function call(problem::Problem, element::Element, ip, time::Float64, ::Type{Val{
     haskey(element, "geometry") || return nothing
     return element("geometry", ip, time)
 end
-

src/problems_heat.jl

@@ -253,4 +253,3 @@ function assemble!{E<:Heat2DSurfaceElements}(assembly::Assembly, problem::Proble
     add!(assembly.K, gdofs, gdofs, K)
     add!(assembly.f, gdofs, fq)
 end
-

src/problems_mortar_2d.jl

@@ -175,7 +175,7 @@ function assemble!(problem::Problem{Mortar}, time::Float64, ::Type{Val{1}}, ::Ty
                 n_s = N1*n1 # normal direction in gauss point
                 xi_m = project_from_slave_to_master(master_element, X_s, n_s, time)
                 N2 = vec(get_basis(master_element, xi_m, time))
-                X_m = N2*X2 
+                X_m = N2*X2
                 De += w*Phi*N1'
                 Me += w*Phi*N2'
                 if props.adjust
@@ -194,7 +194,7 @@ function assemble!(problem::Problem{Mortar}, time::Float64, ::Type{Val{1}}, ::Ty
             # 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]
@@ -210,4 +210,3 @@ function assemble!(problem::Problem{Mortar}, time::Float64, ::Type{Val{1}}, ::Ty
     end # slave elements done, contact virtual work ready
 
 end
-

src/solvers.jl

@@ -10,19 +10,21 @@ type Solver{S<:AbstractSolver}
     norms :: Vector{Tuple}       # solution norms for convergence studies
     ndofs :: Int                 # number of degrees of freedom in problem
     xdmf :: Nullable{Xdmf}       # input/output handle
+    initialized :: Bool
+    u :: Vector{Float64}
+    la :: Vector{Float64}
     properties :: S
 end
 
 
 function Solver{S<:AbstractSolver}(::Type{S}, name="solver", properties...)
     variant = S(properties...)
-    solver = Solver{S}(name, 0.0, [], [], 0, nothing, variant)
+    solver = Solver{S}(name, 0.0, [], [], 0, nothing, false, [], [], variant)
     return solver
 end
 
 function Solver{S<:AbstractSolver}(::Type{S}, problems::Problem...)
-    variant = S()
-    solver = Solver{S}("$(S)Solver", 0.0, [], [], 0, nothing, variant)
+    solver = Solver(S, "$(S)Solver")
     push!(solver.problems, problems...)
     return solver
 end
@@ -237,13 +239,13 @@ Solve linear system using LDLt factorization (SuiteSparse). This version
 requires that final system is symmetric and positive definite, so boundary
 conditions are first eliminated before solution.
 """
-function solve!(K, C1, C2, D, f, g, u, la, ::Type{Val{1}}; F=nothing, debug=false)
+function solve!(solver::Solver, K, C1, C2, D, f, g, u, la, ::Type{Val{1}}; debug=false)
 
-    nnz(D) == 0 || return F, false
+    nnz(D) == 0 || return false
     nz = get_nonzero_rows(C2)
     B = get_nonzero_rows(C2')
     # C2^-1 exists or this doesn't work
-    length(nz) == length(B) || return F, false
+    length(nz) == length(B) || return false
 
     A = get_nonzero_rows(K)
     I = setdiff(A, B)
@@ -270,38 +272,30 @@ function solve!(K, C1, C2, D, f, g, u, la, ::Type{Val{1}}; F=nothing, debug=fals
     end
 
     # solve interior domain using LDLt factorization
-    if F == nothing
-        F = ldltfact(K[I,I])
-    end
+    F = ldltfact(K[I,I])
     u[I] = F \ (f[I] - K[I,B]*u[B])
     # solve lambda
     la[B] = lufact(C1[B,nz]) \ full(f[B] - K[B,I]*u[I] - K[B,B]*u[B])
 
-    return F, true
+    return true
 end
 
 """
 Solve linear system using LU factorization (UMFPACK). This version solves
 directly the saddle point problem without elimination of boundary conditions.
 """
-function solve!(K, C1, C2, D, f, g, u, la, ::Type{Val{2}}; F=nothing)
+function solve!(solver::Solver, K, C1, C2, D, f, g, u, la, ::Type{Val{2}})
     # construct global system Ax = b and solve using lufact (UMFPACK)
     A = [K C1'; C2  D]
     b = [f; g]
-    nz = get_nonzero_rows(A)
-    x = zeros(length(b))
-    if F == nothing
-        F = lufact(A[nz,nz])
-    end
-    x[nz] = F \ full(b[nz])
-    ndofs = size(K, 1)
-    u[:] = x[1:ndofs]
-    la[:] = x[ndofs+1:end]
-    return F, true
+    x = lufact(A) \ full(b)
+    u[:] = x[1:solver.ndofs]
+    la[:] = x[solver.ndofs+1:end]
+    return true
 end
 
 """ Default linear system solver for solver. """
-function solve_linear_system(solver::Solver; F=nothing, empty_assemblies_before_solution=true, show_info=true)
+function solve!(solver::Solver; empty_assemblies_before_solution=true, show_info=true)
     show_info && info("Solving problems ...")
     t0 = Base.time()
 
@@ -317,6 +311,11 @@ function solve_linear_system(solver::Solver; F=nothing, empty_assemblies_before_
     K = 1/2*(K + K')
     M = 1/2*(M + M')
 
+    nz = ones(solver.ndofs)
+    nz[get_nonzero_rows(C2)] = 0.0
+    nz[get_nonzero_rows(D)] = 0.0
+    D += spdiagm(nz)
+
     # free up some memory before solution
     for problem in get_problems(solver)
         if empty_assemblies_before_solution
@@ -327,22 +326,25 @@ function solve_linear_system(solver::Solver; F=nothing, empty_assemblies_before_
         gc()
     end
 
-    u = zeros(solver.ndofs)
-    la = zeros(solver.ndofs)
-
+    ndofs = solver.ndofs
+    u = zeros(ndofs)
+    la = zeros(ndofs)
     status = false
     i = 0
     for i in [1, 2]
-        F, status = solve!(K, C1, C2, D, f, g, u, la, Val{i}; F=F)
+        status = solve!(solver, K, C1, C2, D, f, g, u, la, Val{i})
         status && break
     end
     status || error("Failed to solve linear system!")
-
     t1 = round(Base.time()-t0, 2)
     norms = (norm(u), norm(la))
     show_info && info("Solved problems in $t1 seconds using solver $i. Solution norms = $norms.")
     push!(solver.norms, norms)
-    return F, u, la
+
+    solver.u = u
+    solver.la = la
+
+    return
 end
 
 """ Default assembler for solver. """
@@ -398,7 +400,11 @@ end
 
 """ Default initializer for solver. """
 function initialize!(solver::Solver; show_info=true)
-    show_info && info("Initializing problems ...")
+    if solver.initialized
+        show_info && info("initialize!(): solver already initialized")
+        return
+    end
+    show_info && info("Initializing solver ...")
     problems = get_problems(solver)
     length(problems) != 0 || error("Empty solver, add problems to solver using push!")
     t0 = Base.time()
@@ -419,16 +425,18 @@ function initialize!(solver::Solver; show_info=true)
     info("Total number of nodes in problems: $nnodes")
     maxdof = maximum(nnodes)*field_dim
     info("# of max dof (=size of solution vector) is $maxdof")
-    u = zeros(maxdof)
-    la = zeros(maxdof)
+    solver.u = zeros(maxdof)
+    solver.la = zeros(maxdof)
     # TODO: this could be used to initialize elements too...
+    # TODO: cannot initialize to zero always, construct vector from elements.
     for problem in problems
-        problem.assembly.u = u
-        problem.assembly.la = la
+        problem.assembly.u = zeros(maxdof)
+        problem.assembly.la = zeros(maxdof)
         # initialize(problem, ....)
     end
     t1 = round(Base.time()-t0, 2)
-    show_info && info("Initialized problems in $t1 seconds.")
+    show_info && info("Initialized solver in $t1 seconds.")
+    solver.initialized = true
 end
 
 function get_all_elements(solver::Solver)
@@ -472,7 +480,10 @@ function get_temporal_collection(xdmf::Xdmf)
 end
 
 """ Default update for solver. """
-function update!(solver::Solver, u::Vector, la::Vector; show_info=true)
+function update!{S}(solver::Solver{S}; show_info=true)
+    u = solver.u
+    la = solver.la
+
     show_info && info("Updating problems ...")
     t0 = Base.time()
 
@@ -487,88 +498,105 @@ function update!(solver::Solver, u::Vector, la::Vector; show_info=true)
 
     # if io is attached to solver, update hdf / xml also
     if !isnull(solver.xdmf)
-        xdmf = get(solver.xdmf)
-        temporal_collection = get_temporal_collection(xdmf)
-        frame = new_child(temporal_collection, "Grid")
-        new_child(frame, "Time", Dict("Value" => solver.time))
-
-        # save geometry
-        X = solver("geometry", solver.time)
-        node_ids = sort(collect(keys(X)))
-        geometry = hcat([X[nid] for nid in node_ids]...)
-        ndim, nnodes = size(geometry)
-        geom_type = ndim == 2 ? "XY" : "XYZ"
-        dataitem = new_dataitem(xdmf, "/Node IDs", node_ids)
-        geom = new_child(frame, "Geometry", Dict("Type" => geom_type))
-        dataitem = new_dataitem(xdmf, "/Geometry", geometry)
-        add_child(geom, dataitem)
-
-        # save topology
-        all_elements = get_all_elements(solver)
-        nelements = length(all_elements)
-        element_types = unique(map(get_element_type, all_elements))
-
-        xdmf_element_mapping = Dict(
-            "Seg2" => "Polyline",
-            "Tri3" => "Triangle",
-            "Quad4" => "Quadrilateral",
-            "Tet4" => "Tetrahedron",
-            "Pyramid5" => "Pyramid",
-            "Wedge6" => "Wedge",
-            "Hex8" => "Hexahedron",
-            "Seg3" => "Edge_3",
-            "Tri6" => "Tri_6",
-            "Quad8" => "Quad_8",
-            "Tet10" => "Tet_10",
-            "Pyramid13" => "Pyramid_13",
-            "Wedge15" => "Wedge_15",
-            "Hex20" => "Hex_20")
-
-        for element_type in element_types
-            elements = filter_by_element_type(element_type, all_elements)
-            sort!(elements, by=get_element_id)
-            element_ids = map(get_element_id, elements)
-            element_conn = map(get_connectivity, elements)
-            element_conn = transpose(hcat(element_conn...)) - 1
-            element_code = split(string(element_type), ".")[end]
-            dataitem = new_dataitem(xdmf, "/Topology/$element_code/Element IDs", element_ids)
-            dataitem = new_dataitem(xdmf, "/Topology/$element_code/Connectivity", element_conn)
-            topology = new_child(frame, "Topology")
-            set_attribute(topology, "TopologyType", xdmf_element_mapping[element_code])
-            set_attribute(topology, "NumberOfElements", length(elements))
-            add_child(topology, dataitem)
-        end
-
-        # save solved fields
-        time = "Time $(solver.time)"
-        unknown_field_name = get_unknown_field_name(solver)
-        U = solver(unknown_field_name, solver.time)
-        node_ids2 = sort(collect(keys(U)))
-
-        @assert node_ids == node_ids2
-        ndim = length(U[first(node_ids)])
-        field_type = ndim == 1 ? "Scalar" : "Vector"
-        field_center = "Node"
-        if ndim == 2
-            for nid in node_ids
-                U[nid] = [U[nid]; 0.0]
-            end
-            ndim = 3
-        end
-        U = hcat([U[nid] for nid in node_ids]...)
-        unknown_field_name = ucfirst(unknown_field_name)
-        dataitem = new_dataitem(xdmf, "/Results/$time/Nodal Fields/$unknown_field_name", U)
-        attribute = new_child(frame, "Attribute")
-        set_attribute(attribute, "Name", unknown_field_name)
-        set_attribute(attribute, "Center", field_center)
-        add_child(attribute, dataitem)
-        save!(xdmf)
+        update_xdmf!(solver)
     end
 
     t1 = round(Base.time()-t0, 2)
     show_info && info("Updated problems in $t1 seconds.")
 end
 
+function update_xdmf!{S}(solver::Solver{S}; show_info=true)
+    xdmf = get(solver.xdmf)
+    temporal_collection = get_temporal_collection(xdmf)
+
+    frame = new_element("Grid")
+    new_child(frame, "Time", Dict("Value" => solver.time))
+
+    # save geometry
+    X = solver("geometry", solver.time)
+    node_ids = sort(collect(keys(X)))
+    geometry = hcat([X[nid] for nid in node_ids]...)
+    ndim, nnodes = size(geometry)
+    geom_type = ndim == 2 ? "XY" : "XYZ"
+    dataitem = new_dataitem(xdmf, "/Node IDs", node_ids)
+    geom = new_child(frame, "Geometry", Dict("Type" => geom_type))
+    dataitem = new_dataitem(xdmf, "/Geometry", geometry)
+    add_child(geom, dataitem)
+
+    # save topology
+    all_elements = get_all_elements(solver)
+    nelements = length(all_elements)
+    element_types = unique(map(get_element_type, all_elements))
+
+    xdmf_element_mapping = Dict(
+        "Seg2" => "Polyline",
+        "Tri3" => "Triangle",
+        "Quad4" => "Quadrilateral",
+        "Tet4" => "Tetrahedron",
+        "Pyramid5" => "Pyramid",
+        "Wedge6" => "Wedge",
+        "Hex8" => "Hexahedron",
+        "Seg3" => "Edge_3",
+        "Tri6" => "Tri_6",
+        "Quad8" => "Quad_8",
+        "Tet10" => "Tet_10",
+        "Pyramid13" => "Pyramid_13",
+        "Wedge15" => "Wedge_15",
+        "Hex20" => "Hex_20")
+
+    for element_type in element_types
+        elements = filter_by_element_type(element_type, all_elements)
+        sort!(elements, by=get_element_id)
+        element_ids = map(get_element_id, elements)
+        element_conn = map(get_connectivity, elements)
+        element_conn = transpose(hcat(element_conn...)) - 1
+        element_code = split(string(element_type), ".")[end]
+        dataitem = new_dataitem(xdmf, "/Topology/$element_code/Element IDs", element_ids)
+        dataitem = new_dataitem(xdmf, "/Topology/$element_code/Connectivity", element_conn)
+        topology = new_child(frame, "Topology")
+        set_attribute(topology, "TopologyType", xdmf_element_mapping[element_code])
+        set_attribute(topology, "NumberOfElements", length(elements))
+        add_child(topology, dataitem)
+    end
+
+    # save solved fields
+    unknown_field_name = get_unknown_field_name(solver)
+    U = solver(unknown_field_name, solver.time)
+    node_ids2 = sort(collect(keys(U)))
+
+    @assert node_ids == node_ids2
+    ndim = length(U[first(node_ids)])
+    field_type = ndim == 1 ? "Scalar" : "Vector"
+    field_center = "Node"
+    if ndim == 2
+        for nid in node_ids
+            U[nid] = [U[nid]; 0.0]
+        end
+        ndim = 3
+    end
+    U = hcat([U[nid] for nid in node_ids]...)
+    unknown_field_name = ucfirst(unknown_field_name)
+    time = solver.time
+    path = ""
+    if S == Nonlinear
+        iteration = solver.properties.iteration
+        path = "/Results/Time $time/Iteration $iteration/Nodal Fields/$unknown_field_name"
+    elseif S == Linear
+        path = "/Results/Time $time/Nodal Fields/$unknown_field_name"
+    end
+    dataitem = new_dataitem(xdmf, path, U)
+    attribute = new_child(frame, "Attribute")
+    set_attribute(attribute, "Name", unknown_field_name)
+    set_attribute(attribute, "Center", field_center)
+    set_attribute(attribute, "AttributeType", field_type)
+    add_child(attribute, dataitem)
+    if (S == Linear) || ((S == Nonlinear) && has_converged(solver))
+        add_child(temporal_collection, frame)
+    end
+    save!(xdmf)
+end
+
+
 ### Nonlinear quasistatic solver
 
 type Nonlinear <: AbstractSolver
@@ -580,7 +608,7 @@ type Nonlinear <: AbstractSolver
 end
 
 function Nonlinear()
-    solver = Nonlinear(0, 1, 20, 5.0e-5, true)
+    solver = Nonlinear(0, 1, 10, 5.0e-5, true)
     return solver
 end
 
@@ -646,12 +674,10 @@ function call(solver::Solver{Nonlinear}; show_info=true)
 
         # 2.1 update linearized assemblies
         assemble!(solver)
-
         # 2.2 call solver for linearized system
-        F, u, la = solve_linear_system(solver)
-
+        solve!(solver)
         # 2.3 update solution back to elements
-        update!(solver, u, la)
+        update!(solver)
 
         # 2.4 check convergence
         if has_converged(solver)
@@ -713,7 +739,7 @@ function assemble!(solver::Solver{Linear}; show_info=true)
     show_info && info("Assembled $nproblems problems in $t1 seconds. ndofs = $ndofs.")
 end
 
-function call(solver::Solver{Linear}; F=nothing, show_info=true, return_factorization=true)
+function call(solver::Solver{Linear}; show_info=true)
     t0 = Base.time()
     show_info && info(repeat("-", 80))
     show_info && info("Starting linear solver")
@@ -721,13 +747,10 @@ function call(solver::Solver{Linear}; F=nothing, show_info=true, return_factoriz
     show_info && info(repeat("-", 80))
     initialize!(solver)
     assemble!(solver)
-    F, u, la = solve_linear_system(solver; F=F, empty_assemblies_before_solution=false)
-    update!(solver, u, la)
+    solve!(solver)
+    update!(solver)
     t1 = round(Base.time()-t0, 2)
     show_info && info("Linear solver ready in $t1 seconds.")
-    if return_factorization
-        return F
-    end
 end
 
 """ Convenience function to call linear solver. """