Add postprocessing features (#100)

* refactored code for solvers.

* Added elementary tests for least-squares fitting of strain and stress fields

* A more realistic postprocess + Xdmf writing test

* removed debug keyword argument from test

* Rewrite update_xdmf!

New function to update Xdmf file no longer takes Solver object but
xdmf, problem, time and fields to write, for example

julia> update_xdmf!(xdmf, problem, 0.0, ["displacement", "temperature"])

All problems are written separately and put together into one
SpatialCollection, allowing to have more structured Xdmf and making
it easier to write complicated field configurations. Support for Xdmf
API 3.0 added.

* Support for Tensor6 field writing

* moved update_xdmf! to io.jl

* Removed some empty files

* Not use old Postprocessor, obsolete code.

* Not use old XDMF (obsolete code). Fixed test.

* removed some postprocessing to pass test, maybe we should drop abaqus.jl from code as obsolete

* add function get_temporal_collection back, it's used by update_xdmf of modal solver

* postprocess of boundary problems also

* added test for contact pressure. dl+quad test output was written in wrong file, fixed.

* postprocess for contact pressure

* contact pressure postprocess

* with boundary problems always store also the primary unknown field

* Change "reaction force" -> "lambda"

* testing postprocess of reaction force also

* sign convention

src/JuliaFEM.jl

@@ -30,9 +30,6 @@ module Testing
 
 end
 
-include("io.jl")
-export Xdmf, h5file, xmffile, xdmf_filter, new_dataitem, get_temporal_collection
-
 include("fields.jl")
 export Field, DCTI, DVTI, DCTV, DVTV, CCTI, CVTI, CCTV, CVTV, Increment
 
@@ -42,7 +39,8 @@ export AbstractPoint, Point, IntegrationPoint, IP, Node
 ### ELEMENTS ###
 include("elements.jl") # common element routines
 export Node, AbstractElement, Element, update!, get_connectivity, get_basis,
-       get_dbasis, inside, get_local_coordinates
+       get_dbasis, inside, get_local_coordinates, get_element_type,
+       filter_by_element_type, get_element_id
 
 include("elements_lagrange.jl") # Continuous Galerkin (Lagrange) elements
 export get_reference_coordinates,
@@ -99,6 +97,8 @@ export calculate_normals, calculate_normals!, project_from_slave_to_master,
        project_from_master_to_slave, Mortar, get_slave_elements,
        get_polygon_clip
 
+include("io.jl")
+export Xdmf, h5file, xmffile, xdmf_filter, new_dataitem, update_xdmf!, save!
 
 ### ASSEMBLY + SOLVE ###
 include("assembly.jl")
@@ -141,17 +141,11 @@ export aster_create_elements, parse_aster_med_file, is_aster_mail_keyword,
 end
 
 module Postprocess
-
 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
-
-include("postprocess_xdmf.jl")
-export XDMF, xdmf_new_result!, xdmf_save_field!, xdmf_save!, DataFrame
-
 end
-export Postprocessor
 
 module Abaqus
 include("abaqus.jl")

src/abaqus.jl

@@ -564,8 +564,8 @@ function process_output_request(model::Model, solver::Solver, output_request::Ab
             haskey(code_mapping, code) || continue
             field_name = code_mapping[code]
             abbr = get(abbr_mapping, code, code)
-            table = solver(DataFrame, field_name, abbr, solver.time)
-            push!(tables, table)
+            #table = solver(DataFrame, field_name, abbr, solver.time)
+            #push!(tables, table)
         end
         length(tables) != 0 || continue
         results = join(tables..., on=:NODE, kind=:outer)

src/api.jl

@@ -1,5 +0,0 @@
-# This file is a part of JuliaFEM.
-# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
-
-include("api_types.jl")
-include("api_functions.jl")

src/io.jl

@@ -8,6 +8,8 @@ type Xdmf
     name :: String
     xml :: XMLElement
     hdf :: HDF5File
+    hdf_counter :: Int
+    format :: String
 end
 
 function Xdmf()
@@ -23,7 +25,7 @@ function xmffile(xdmf::Xdmf)
 end
 
 """ Initialize a new Xdmf object. """
-function Xdmf(name::String; overwrite=false)
+function Xdmf(name::String; version="3.0", overwrite=false)
     xdmf = new_element("Xdmf")
     h5file = "$name.h5"
     xmlfile = "$name.xmf"
@@ -47,16 +49,14 @@ function Xdmf(name::String; overwrite=false)
     end
 
     set_attribute(xdmf, "xmlns:xi", "http://www.w3.org/2001/XInclude")
-    set_attribute(xdmf, "Version", "2.1")
+    set_attribute(xdmf, "Version", version)
     flag = isfile(h5file) ? "r+" : "w"
     hdf = h5open(h5file, flag)
-    return Xdmf(name, xdmf, hdf)
+    return Xdmf(name, xdmf, hdf, 1, "HDF")
 end
 
 """ Return the basic structure of Xdmf document. Creates a new TemporalCollection if not found.
-
 Basic structure for XML part of Xdmf file is
-
 <?xml version="1.0" encoding="utf-8"?>
 <Xdmf xmlns:xi="http://www.w3.org/2001/XInclude" Version="2.1">
   <Domain>
@@ -64,7 +64,6 @@ Basic structure for XML part of Xdmf file is
     </Grid>
   </Domain>
 </Xdmf>
-
 """
 function get_temporal_collection(xdmf::Xdmf)
     domain = find_element(xdmf.xml, "Domain")
@@ -194,6 +193,13 @@ function traverse(xdmf::Xdmf, x::XMLElement, attr_name::String)
     if '/' in attr_name
         items = split(attr_name, '/')
         new_item = xdmf_filter(childs, first(items))
+        if new_item == nothing
+            info("traverse: childs:")
+            for child in childs
+                info(LightXML.name(child))
+            end
+            error("traverse: failed, items = $items, xdmf_filter not find child")
+        end
         new_path = join(items[2:end], '/')
         return traverse(xdmf, new_item, new_path)
     end
@@ -211,11 +217,14 @@ function read(xdmf::Xdmf, path::String)
     result = traverse(xdmf, xdmf.xml, path)
     if endswith(path, "DataItem")
         format = attribute(result, "Format"; required=true)
-        @assert format == "HDF"
-        h5file, path = map(String, split(content(result), ':'))
-        h5file = dirname(xdmf.name) * "/" * h5file
-        isfile(h5file) || throw("Xdmf: h5 file $h5file not found!")
-        return read(xdmf.hdf, path)
+        if format == "HDF"
+            h5file, path = map(String, split(content(result), ':'))
+            h5file = dirname(xdmf.name) * "/" * h5file
+            isfile(h5file) || throw("Xdmf: h5 file $h5file not found!")
+            return read(xdmf.hdf, path)
+        else
+            error("Read from Xdmf, reading from $format not implemented")
+        end
     else
         return result
     end
@@ -228,14 +237,14 @@ function save!(xdmf::Xdmf)
     save_file(doc, xmffile(xdmf))
 end
 
-function new_dataitem{T,N}(xdmf::Xdmf, path::String, data::Array{T,N}; format="HDF")
+function new_dataitem{T,N}(xdmf::Xdmf, path::String, data::Array{T,N})
     dataitem = new_element("DataItem")
     datatype = replace("$T", "64", "")
     dimensions = join(reverse(size(data)), " ")
     set_attribute(dataitem, "DataType", datatype)
     set_attribute(dataitem, "Dimensions", dimensions)
-    set_attribute(dataitem, "Format", format)
-    if format == "HDF"
+    set_attribute(dataitem, "Format", xdmf.format)
+    if xdmf.format == "HDF"
         hdf = basename(h5file(xdmf))
         if exists(xdmf.hdf, path)
             info("Xdmf: $path already existing in h5 file, not overwriting.")
@@ -243,9 +252,194 @@ function new_dataitem{T,N}(xdmf::Xdmf, path::String, data::Array{T,N}; format="H
             write(xdmf.hdf, path, data)
         end
         add_text(dataitem, "$hdf:$path")
-    elseif format == "XML"
-        add_text(dataitem, strip(string(data), ['[', ']']))
+    elseif xdmf.format == "XML"
+        text_data = string(data')
+        text_data = strip(text_data, ['[', ']'])
+        text_data = replace(text_data, ';', '\n')
+        text_data = "\n" * text_data * "\n"
+        add_text(dataitem, text_data)
+    else
+        error("Unsupported Xdmf big data format $(xdmf.format)")
     end
     return dataitem
 end
 
+""" Create a new DataItem element, hdf path automatically determined. """
+function new_dataitem{T,N}(xdmf::Xdmf, data::Array{T,N})
+    if xdmf.format == "XML"
+        # Path can be whatever as XML format does not store to HDF at all
+        return new_dataitem(xdmf, "/whatever", data)
+    else
+        debug("Determining path for HDF file automatically.")
+        path = "/DataItem_$(xdmf.hdf_counter)"
+        while exists(xdmf.hdf, path)
+            xdmf.hdf_counter += 1
+            path = "/DataItem_$(xdmf.hdf_counter)"
+        end
+        debug("HDF path automatically determined to be $path")
+        return new_dataitem(xdmf, path, data)
+    end
+end
+
+global const xdmf_element_mapping = Dict(
+        "Poi1" => "Polyvertex",
+        "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")
+
+""" Write new fields to Xdmf file.
+
+Examples
+--------
+
+To write displacement and temperature fields from p1 at time t=0.0:
+
+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")
+
+    # 1. find domain
+    xml = xdmf.xml
+    domain = find_element(xml, "Domain")
+    if domain == nothing
+        info("Xdmf: Domain not found, creating.")
+        domain = new_child(xml, "Domain")
+    else
+        debug("Xdmf: Domain already defined, skipping.")
+    end
+
+    # 2. find for TemporalCollection
+    temporal_collection = find_element(domain, "Grid")
+    if temporal_collection == nothing
+        info("Xdmf: Temporal collection not found, creating.")
+        temporal_collection = new_child(domain, "Grid")
+        set_attribute(temporal_collection, "GridType", "Collection")
+        set_attribute(temporal_collection, "Name", "Time")
+        set_attribute(temporal_collection, "CollectionType", "Temporal")
+    else
+        debug("Xdmf: Temporal collection found, skipping.")
+    end
+
+    # 2.1 make sure that Grid element we found really is TemporalCollection
+    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"
+
+    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.")
+            return
+        end
+    end
+
+    frame_name = problem.name
+    info("Xdmf: Creating Grid for problem $frame_name")
+    frame = new_child(spatial_collection, "Grid")
+    set_attribute(frame, "Name", frame_name)
+
+    # 4. save geometry
+    X_dict = problem("geometry", time)
+    node_ids = sort(collect(keys(X_dict)))
+    node_mapping = Dict(j => i for (i, j) in enumerate(node_ids))
+    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")
+    X_dataitem = new_dataitem(xdmf, X_array)
+    geometry = new_child(frame, "Geometry")
+    set_attribute(geometry, "Type", geom_type)
+    add_child(geometry, X_dataitem)
+
+    # 5. save topology
+    all_elements = get_elements(problem)
+    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.")
+
+    for element_type in element_types
+        elements = filter_by_element_type(element_type, all_elements)
+        nelements = length(elements)
+        info("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)
+        element_conn = hcat(element_conn...)
+        element_code = split(string(element_type), ".")[end]
+        topology_dataitem = new_dataitem(xdmf, 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, topology_dataitem)
+    end
+
+    # 6. save requested fields
+    for field_name in fields
+        field_dict = problem(field_name, time)
+        field_center = "Node"
+        field_node_ids = sort(collect(keys(field_dict)))
+        @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")
+            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")
+
+        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, "Center", field_center)
+        set_attribute(attribute, "AttributeType", field_type)
+        add_child(attribute, field_dataitem)
+    end
+
+    save!(xdmf)
+    info("Xdmf: all done.")
+end

src/postprocess_xdmf.jl

@@ -1,194 +0,0 @@
-# This file is a part of JuliaFEM.
-# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
-
-using JuliaFEM
-using LightXML
-
-# element codes: http://www.paraview.org/pipermail/paraview/2013-July/028859.html
-# > from  ./VTK/ThirdParty/xdmf2/vtkxdmf2/libsrc/XdmfTopology.h
-# >
-# > // Topologies
-# > #define XDMF_NOTOPOLOGY     0x0
-# > #define XDMF_POLYVERTEX     0x1
-# > #define XDMF_POLYLINE       0x2
-# > #define XDMF_POLYGON        0x3
-# > #define XDMF_TRI            0x4
-# > #define XDMF_QUAD           0x5
-# > #define XDMF_TET            0x6
-# > #define XDMF_PYRAMID        0x7
-# > #define XDMF_WEDGE          0x8
-# > #define XDMF_HEX            0x9
-# > #define XDMF_EDGE_3         0x0022
-# > #define XDMF_TRI_6          0x0024
-# > #define XDMF_QUAD_8         0x0025
-# > #define XDMF_QUAD_9         0x0023
-# > #define XDMF_TET_10         0x0026
-# > #define XDMF_PYRAMID_13     0x0027
-# > #define XDMF_WEDGE_15       0x0028
-# > #define XDMF_WEDGE_18       0x0029
-# > #define XDMF_HEX_20         0x0030
-# > #define XDMF_HEX_24         0x0031
-# > #define XDMF_HEX_27         0x0032
-# > #define XDMF_MIXED          0x0070
-# > #define XDMF_2DSMESH        0x0100
-# > #define XDMF_2DRECTMESH     0x0101
-# > #define XDMF_2DCORECTMESH   0x0102
-# > #define XDMF_3DSMESH        0x1100
-# > #define XDMF_3DRECTMESH     0x1101
-# > #define XDMF_3DCORECTMESH   0x1102
-
-get_xdmf_element_code(element::Element{Poi1})  = 0x0001
-get_xdmf_element_code(element::Element{Seg2})  = 0x0002
-get_xdmf_element_code(element::Element{Seg3})  = 0x0003
-get_xdmf_element_code(element::Element{Tri3})  = 0x0004
-get_xdmf_element_code(element::Element{Quad4}) = 0x0005
-get_xdmf_element_code(element::Element{Tet4})  = 0x0006
-get_xdmf_element_code(element::Element{Hex8})  = 0x0009
-get_xdmf_element_code(element::Element{Tet10}) = 0x0026
-
-type XDMF
-    dimension :: Int
-    use_hdf :: Bool
-    xdoc :: XMLDocument
-    domain :: XMLElement
-    temporal_collection :: XMLElement
-    current_grid
-    permutation :: Vector{Int}
-end
-
-function XDMF()
-    xdoc = XMLDocument()
-    xroot = create_root(xdoc, "Xdmf")
-    set_attribute(xroot, "xmlns:xi", "http://www.w3.org/2001/XInclude")
-    set_attribute(xroot, "Version", "2.1")
-    domain = new_child(xroot, "Domain")
-    temporal_collection = new_child(domain, "Grid")
-    set_attribute(temporal_collection, "CollectionType", "Temporal")
-    set_attribute(temporal_collection, "GridType", "Collection")
-    set_attribute(temporal_collection, "Name", "Collection")
-    return XDMF(3, false, xdoc, domain, temporal_collection, Union{}, [])
-end
-
-function xdmf_new_result!(xdmf::XDMF, elements::Vector, time)
-    grid = new_child(xdmf.temporal_collection, "Grid")
-    set_attribute(grid, "Name", "Grid")
-    time_ = new_child(grid, "Time")
-    set_attribute(time_, "Value", time)
-    xdmf.current_grid = grid
-
-    # 1. calculate permutation
-    nids = Set()
-    X = Dict{Int64, Vector{Float64}}()
-    for element in elements
-        conn = get_connectivity(element)
-        push!(nids, conn...)
-        X_el = element["geometry"](time)
-        for (i, c) in enumerate(conn)
-            X[c] = X_el[i]
-        end
-    end
-    xdmf.permutation = sort(collect(nids))
-    iperm = Dict{Int64, Int64}()
-    for (i, j) in enumerate(xdmf.permutation)
-        iperm[j] = i
-    end
-
-    # 2. write nodes
-    geometry = new_child(grid, "Geometry")
-    set_attribute(geometry, "Type", xdmf.dimension == 3 ? "XYZ" : "XY")
-    dataitem = new_child(geometry, "DataItem")
-    set_attribute(dataitem, "DataType", "Float")
-    set_attribute(dataitem, "Format", "XML")
-    #set_attribute(dataitem, "Precision", 8)
-    s = []
-    ndim = 0
-    for i in xdmf.permutation
-        ndim += length(X[i])
-        push!(s, join(round(X[i], 5), " "))
-    end
-    set_attribute(dataitem, "Dimensions", ndim)
-    add_text(dataitem, "\n"*join(s, "\n")*"\n")
-
-    # 3. write elements
-    topology = new_child(grid, "Topology")
-    set_attribute(topology, "TopologyType", "Mixed")
-    set_attribute(topology, "NumberOfElements", length(elements))
-    dataitem = new_child(topology, "DataItem")
-    set_attribute(dataitem, "Format", "XML")
-    set_attribute(dataitem, "DataType", "Int")
-#   set_attribute(dataitem, "Precision", 8)
-    s = []
-    eldim = 0
-    for element in elements
-        eltype = get_xdmf_element_code(element)
-        # note: id numbers start from 0 in Xdmf
-        conn = [iperm[j] for j in get_connectivity(element)] - 1
-        data = [eltype; conn]
-        eldim += length(data)
-        push!(s, join(data, " "))
-    end
-    set_attribute(dataitem, "Dimensions", eldim)
-    add_text(dataitem, "\n"*join(s, "\n")*"\n")
-end
-
-function xdmf_save_field!(xdmf, elements::Vector, time, field_name; field_type="Scalar", debug=false)
-    f = Dict()
-    field_dim = 0
-    for element in elements
-        haskey(element, field_name) || continue
-        g = element[field_name](time)
-        conn = get_connectivity(element)
-        for (i, c) in enumerate(conn)
-            gi = g[i]
-            if (field_type == "Vector") && (length(gi) < 3)
-                # paraview goes crazy if 2d model with 2d displacement vector
-                gi = [gi; 0.0]
-            end
-            if field_dim == 0
-                field_dim = length(gi)
-            end
-            field_dim == length(gi) || error("several dimensions in field, dim = $field_dim.")
-            f[c] = gi
-        end
-    end
-
-    if length(f) == 0
-        warn("xdmf_save_field!(): field $field_name was not found from set of elements")
-        return
-    end
-
-    attribute = new_child(xdmf.current_grid, "Attribute")
-    set_attribute(attribute, "Center", "Node")
-    set_attribute(attribute, "Name", ucfirst(field_name))
-    set_attribute(attribute, "Type", field_type)
-    dataitem = new_child(attribute, "DataItem")
-    set_attribute(dataitem, "DataType", "Float")
-    set_attribute(dataitem, "Format", "XML")
-    #set_attribute(dataitem, "Precision", 8)
-    debug && info("field dim = $field_dim")
-    debug && info(f)
-    s = []
-    dim = 0
-    for i in xdmf.permutation
-        gi = zeros(field_dim)
-        if haskey(f, i)
-            gi = f[i]
-        end
-        push!(s, join(round(gi, 5), " "))
-        dim += length(gi)
-    end
-    set_attribute(dataitem, "Dimensions", dim)
-    add_text(dataitem, "\n"*join(s, "\n")*"\n")
-end
-
-function xdmf_save_field!(xdmf, problem::Problem, time, field_name; field_type="Scalar")
-    xdmf_save_field!(xdmf, problem.elements, time, field_name; field_type=field_type)
-end
-
-function xdmf_new_result!(xdmf, problem::Problem, time)
-    xdmf_new_result!(xdmf, problem.elements, time)
-end
-
-function xdmf_save!(xdmf, filename)
-    save_file(xdmf.xdoc, filename)
-end

src/problems.jl

@@ -89,6 +89,7 @@ type Problem{P<:AbstractProblem}
     dofmap :: Dict{Element, Vector{Int64}} # connects element local dofs to global dofs
     assembly :: Assembly
     fields :: Dict{AbstractString, Field}
+    postprocess_fields :: Vector{String}
     properties :: P
 end
 
@@ -103,10 +104,10 @@ julia> prob2 = Problem(Elasticity, 3)
 
 """
 function Problem{P<:FieldProblem}(::Type{P}, name::AbstractString, dimension::Int64)
-    return Problem{P}(name, dimension, "none", [], Dict(), Assembly(), Dict(), P())
+    return Problem{P}(name, dimension, "none", [], Dict(), Assembly(), Dict(), Vector(), P())
 end
 function Problem{P<:FieldProblem}(::Type{P}, dimension::Int64)
-    return Problem{P}("$P problem", dimension, "none", [], Dict(), Assembly(), Dict(), P())
+    return Problem(P, "$P problem", dimension)
 end
 
 """ Construct a new boundary problem.
@@ -119,28 +120,27 @@ 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(), Dict(), P())
+    return Problem{P}(name, dimension, parent_field_name, [], Dict(), Assembly(), Dict(), Vector(), 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(), Dict(), P())
+    return Problem{P}(name, dimension, parent_field_name, [], Dict(), Assembly(), Dict(), Vector(), P())
 end
 
-function get_formulation_type{P<:FieldProblem}(problem::Problem{P})
+function get_formulation_type(problem::Problem)
     return :incremental
 end
 
-function get_formulation_type{P<:BoundaryProblem}(problem::Problem{P})
-    return :incremental
+function get_unknown_field_name{P<:BoundaryProblem}(::Type{P})
+    return "lambda"
 end
 
 function get_assembly(problem)
     return problem.assembly
 end
 
-
 """ Initialize element ready for calculation. """
 function initialize!(problem::Problem, element::Element, time::Float64)
     field_name = get_unknown_field_name(problem)
@@ -175,7 +175,7 @@ function initialize!(problem::Problem, time::Float64=0.0)
 end
 
 """ Update problem solution vector for assembly. """
-function update!(problem::Problem, assembly::Assembly, u::Vector, la::Vector; verbose=false)
+function update!(problem::Problem, assembly::Assembly, u::Vector, la::Vector)
 
     # resize & fill with zeros vectors if length mismatch with current solution
 
@@ -186,7 +186,7 @@ function update!(problem::Problem, assembly::Assembly, u::Vector, la::Vector; ve
     end
 
     if length(la) != length(assembly.la)
-        info("resizing lagrange multipliers vector u")
+        info("resizing lagrange multiplier vector la")
         resize!(assembly.la, length(la))
         fill!(assembly.la, 0.0)
     end
@@ -199,15 +199,12 @@ function update!(problem::Problem, assembly::Assembly, u::Vector, la::Vector; ve
     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
         assembly.la = la
     elseif get_formulation_type(problem) == :incremental
-        verbose && info("$(problem.name): incremental formulation, adding increment to solution vector")
         assembly.u += u
         assembly.la = la
     elseif get_formulation_type(problem) == :forwarddiff
-        verbose && info("$(problem.name): forwarddiff formulation, adding increment to solution vector and reaction force vector")
         assembly.u += u
         assembly.la += la
     else
@@ -259,13 +256,12 @@ end
 function update!{P<:BoundaryProblem}(problem::Problem{P}, assembly::Assembly, elements::Vector{Element}, time::Float64)
     u, la = get_global_solution(problem, assembly)
     parent_field_name = get_parent_field_name(problem) # displacement
-    field_name = get_unknown_field_name(problem) # reaction force
-    # update solution u and reaction force λ for boundary elements
+    field_name = get_unknown_field_name(problem) # lambda
+    # update solution and lagrange multipliers for boundary elements
     for element in elements
         connectivity = get_connectivity(element)
         update!(element, parent_field_name, time => u[connectivity])
-        # FIXME
-        update!(element, field_name, time => -la[connectivity])
+        update!(element, field_name, time => la[connectivity])
     end
 end
 

src/problems_contact.jl

@@ -58,22 +58,6 @@ function Contact()
         default_fields)
 end
 
-function get_unknown_field_name(problem::Problem{Contact})
-    return "reaction force"
-end
-
-function get_formulation_type(problem::Problem{Contact})
-    #=
-    if problem.properties.use_forwarddiff
-        return :forwarddiff
-    else
-        return :incremental
-    end
-    =#
-    return :incremental
-    #return :forwarddiff
-end
-
 function assemble!(problem::Problem{Contact}, time::Real)
     if problem.properties.dimension == -1
         problem.properties.dimension = dim = size(first(problem.elements), 1)

src/problems_contact_2d.jl

@@ -77,7 +77,7 @@ function assemble!(problem::Problem{Contact}, time::Float64, ::Type{Val{1}}, ::T
         nsl = length(slave_element)
         X1 = slave_element("geometry", time)
         u1 = slave_element("displacement", time)
-        la1 = slave_element("reaction force", time)
+        la1 = slave_element("lambda", time)
         n1 = slave_element("normal", time)
         t1 = slave_element("tangent", time)
         x1 = X1 + u1
@@ -160,8 +160,8 @@ function assemble!(problem::Problem{Contact}, time::Float64, ::Type{Val{1}}, ::T
                 Te += w*reshape(kron(N2, n_s, Phi), 2, 4)
                 He += w*reshape(kron(N1, t_s, Phi), 2, 4)
                 ge += w*Phi*dot(n_s, x_m-x_s)
-                ce += w*N1*dot(n_s, -la_s)
-                Rn += w*dot(n_s, -la_s)
+                ce += w*N1*dot(n_s, la_s)
+                Rn += w*dot(n_s, la_s)
 
                 contact_area += w
                 contact_error += 1/2*w*dot(n_s, x_s-x_m)^2

src/problems_contact_3d.jl

@@ -117,7 +117,7 @@ function assemble!(problem::Problem{Contact}, slave_element::Element{Tri3}, time
     u1 = slave_element("displacement", time)
     x1 = X1 + u1
     n1 = slave_element("normal", time)
-    la = slave_element("reaction force", time)
+    la = slave_element("lambda", time)
 
     Q3 = create_rotation_matrix(slave_element, time)
 
@@ -277,7 +277,7 @@ function assemble!(problem::Problem{Contact}, slave_element::Element{Tri6}, time
             #u1 = sub_slave_element("displacement", time)
             #x1 = X1 + u1
             n1 = sub_slave_element("normal", time)
-            #la = sub_slave_element("reaction force", time)
+            #la = sub_slave_element("lambda", time)
 
             # create auxiliary plane
             xi = mean(get_reference_coordinates(sub_slave_element))
@@ -599,16 +599,16 @@ function assemble!(problem::Problem{Contact}, time::Float64, ::Type{Val{2}}, ::T
         end
         complementarity_condition[j] = contact_pressure[j] - weighted_gap[j]
         
-        if complementarity_condition[j][1] < 0.0
-            is_inactive[j] = 1
-            is_active[j] = 0
-            is_slip[j] = 0
-            is_stick[j] = 0
-        else
+        if complementarity_condition[j][1] > 0.0
             is_inactive[j] = 0
             is_active[j] = 1
             is_slip[j] = 1
             is_stick[j] = 0
+        else
+            is_inactive[j] = 1
+            is_active[j] = 0
+            is_slip[j] = 0
+            is_stick[j] = 0
         end
     end
 
@@ -670,3 +670,18 @@ function assemble!(problem::Problem{Contact}, time::Float64, ::Type{Val{2}}, ::T
     problem.assembly.g = g
 
 end
+
+function postprocess!(problem::Problem{Contact}, time::Float64, ::Type{Val{Symbol("contact pressure")}})
+    n = problem("normal", time)
+    la = problem("lambda", time)
+    node_ids = keys(n)
+    cp = Dict(nid => dot(n[nid], la[nid]) for nid in node_ids)
+    # FIXME: have to define zero contact pressure & lambda to master elements
+    # elements because interface.elements = [slave_elements; master_elements]
+    for nid in keys(la)
+        if !haskey(cp, nid)
+            cp[nid] = 0.0
+        end
+    end
+    update!(problem, "contact pressure", time => cp)
+end

src/problems_dirichlet.jl

@@ -15,10 +15,6 @@ function Dirichlet()
     Dirichlet(:incremental, false, false, 1)
 end
 
-function get_unknown_field_name(::Type{Dirichlet})
-    return "reaction force"
-end
-
 function get_formulation_type(problem::Problem{Dirichlet})
     return problem.properties.formulation
 end
@@ -142,53 +138,8 @@ function assemble!(assembly::Assembly, problem::Problem{Dirichlet},
 
 end
 
-#=
-
-function assemble!(assembly::Assembly, problem::Problem{DirichletProblem}, element::Element, time::Real)
-
-    # get dimension and name of PARENT field
-    field_dim = problem.parent_field_dim
-    field_name = problem.parent_field_name
-
-    gdofs = get_gdofs(element, field_dim)
-    for ip in get_integration_points(element, Val{2})
-        w = ip.weight
-        J = get_jacobian(element, ip, time)
-        JT = transpose(J)
-        if size(JT, 2) == 1  # plane problem
-            w *= norm(JT)
-        else
-            w *= norm(cross(JT[:,1], JT[:,2]))
-        end
-        N = element(ip, time)
-        A = w*N'*N
-
-        if haskey(element, field_name)
-            # add all dimensions at once if defined
-            # element["blaa"] = 0.0
-            # or
-            # element["blaa"] = Vector{Float64}[[0.1, 0.2], [0.3, 0.4]]
-            g = element(field_name, ip, time)
-            if length(g) != length(N)
-                g = g*ones(length(N))
-            end
-            for i=1:field_dim
-                ldofs = gdofs[i:field_dim:end]
-                add!(assembly.C1, ldofs, ldofs, A)
-                add!(assembly.C2, ldofs, ldofs, A)
-            end
-            add!(assembly.g, gdofs, w*g*N)
-        end
-
-        for i=1:field_dim
-            if haskey(element, field_name*" $i")
-                g = element(field_name*" $i", ip, time)
-                ldofs = gdofs[i:field_dim:end]
-                add!(assembly.C1, ldofs, ldofs, A)
-                add!(assembly.C2, ldofs, ldofs, A)
-                add!(assembly.g, ldofs, w*g*N)
-            end
-        end
-    end
+function postprocess!(problem::Problem{Dirichlet}, time::Float64, ::Type{Val{Symbol("reaction force")}})
+    la = problem("lambda", time)
+    rf = Dict(nid => -lai for (nid, lai) in la)
+    update!(problem, "reaction force", time => rf)
 end
-=#

src/problems_elasticity.jl

@@ -534,3 +534,80 @@ function assemble{El<:Elasticity3DSurfaceElements}(problem::Problem{Elasticity},
     end
     return Km, Kg, f
 end
+
+""" Return strain tensor. """
+function get_strain_tensor(problem, element, ip, time)
+    gradu = element("displacement", ip, time, Val{:Grad})
+    eps = 0.5*(gradu' + gradu)
+    return eps
+end
+
+""" Return stress tensor. """
+function get_stress_tensor(problem, element, ip, time)
+    eps = get_strain_tensor(problem, element, ip, time)
+    E = element("youngs modulus", 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
+    return S
+end
+
+""" Return stain vector in "ABAQUS" order 11, 22, 33, 12, 23, 13. """
+function get_strain_vector(problem, element, ip, time)
+    eps = get_strain_tensor(problem, element, ip, time)
+    return [eps[1,1], eps[2,2], eps[3,3], eps[1,2], eps[2,3], eps[1,3]]
+end
+
+""" Return stress vector in "ABAQUS" order 11, 22, 33, 12, 23, 13. """
+function get_stress_vector(problem, element, ip, time)
+    S = get_stress_tensor(problem, element, ip, time)
+    return [S[1,1], S[2,2], S[3,3], S[1,2], S[2,3], S[1,3]]
+end
+
+""" Make least squares fit for some field to nodes. """
+function lsq_fit(problem, elements, field, time)
+    A = SparseMatrixCOO()
+    b = SparseMatrixCOO()
+    volume = 0.0
+    for element in elements
+        gdofs = get_connectivity(element)
+        for ip in get_integration_points(element)
+            detJ = element(ip, time, Val{:detJ})
+            w = ip.weight*detJ
+            N = element(ip, time)
+            f = field(problem, element, ip, time)
+            add!(A, gdofs, gdofs, w*kron(N', N))
+            for i=1:length(f)
+                add!(b, gdofs, w*f[i]*N, i)
+            end
+            volume += w
+        end
+    end
+    debug("Mass matrix for least-squares fit is assembled. Total volume to fit: $volume")
+    A = sparse(A)
+    b = sparse(b)
+    A = 1/2*(A + A')
+    
+    nz = get_nonzero_rows(A)
+    F = ldltfact(A[nz,nz])
+
+    x = F \ b[nz, :]
+
+    nodal_values = Dict(node_id => vec(full(x[idx,:])) for (idx, node_id) in enumerate(nz))
+    return nodal_values
+end
+
+
+""" Postprocessing, extrapolate strain to nodes using least-squares fit. """
+function postprocess!(problem::Problem{Elasticity}, time::Float64, ::Type{Val{:strain}})
+    elements = get_elements(problem)
+    strain = lsq_fit(problem, elements, get_strain_vector, time)
+    update!(elements, "strain", time => strain)
+end
+
+function postprocess!(problem::Problem{Elasticity}, time::Float64, ::Type{Val{:stress}})
+    elements = get_elements(problem)
+    stress = lsq_fit(problem, elements, get_stress_vector, time)
+    update!(elements, "stress", time => stress)
+end

src/problems_mortar.jl

@@ -52,18 +52,6 @@ function Mortar()
     return Mortar(-1, false, false, false, false, Inf, true, true, true, 0.0, 1.0e-9, default_fields)
 end
 
-function get_unknown_field_name(problem::Problem{Mortar})
-    return "reaction force"
-end
-
-function get_formulation_type(problem::Problem{Mortar})
-    if problem.properties.use_forwarddiff
-        return :forwarddiff
-    else
-        return :incremental
-    end
-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!")

src/solvers.jl

@@ -77,7 +77,7 @@ If several field problems exists, they are simply summed together, so
 problems must have unique node ids.
 
 """
-function get_field_assembly(solver::Solver; show_info=true)
+function get_field_assembly(solver::Solver)
     problems = get_field_problems(solver)
 
     M = SparseMatrixCOO()
@@ -96,7 +96,7 @@ function get_field_assembly(solver::Solver; show_info=true)
 
     if solver.ndofs == 0
         solver.ndofs = size(K, 1)
-        show_info && info("automatically determined problem dimension, ndofs = $(solver.ndofs)")
+        info("automatically determined problem dimension, ndofs = $(solver.ndofs)")
     end
 
     M = sparse(M, solver.ndofs, solver.ndofs)
@@ -364,8 +364,8 @@ function solve!(solver::Solver; empty_assemblies_before_solution=true, symmetric
 end
 
 """ Default assembler for solver. """
-function assemble!(solver::Solver; show_info=true, timing=true, with_mass_matrix=false)
-    show_info && info("Assembling problems ...")
+function assemble!(solver::Solver; timing=true, with_mass_matrix=false)
+    info("Assembling problems ...")
 
     function do_assemble(problem)
         t00 = Base.time()
@@ -391,7 +391,7 @@ function assemble!(solver::Solver; show_info=true, timing=true, with_mass_matrix
 
     solver.ndofs = ndofs
     t1 = round(Base.time()-t0, 2)
-    show_info && info("Assembled $nproblems problems in $t1 seconds. ndofs = $ndofs.")
+    info("Assembled $nproblems problems in $t1 seconds. ndofs = $ndofs.")
     if timing
         info("Assembly times:")
         for (i, problem) in enumerate(solver.problems)
@@ -423,12 +423,12 @@ function get_unknown_field_dimension(solver::Solver)
 end
 
 """ Default initializer for solver. """
-function initialize!(solver::Solver; show_info=true)
+function initialize!(solver::Solver)
     if solver.initialized
-        show_info && info("initialize!(): solver already initialized")
+        warn("initialize!(): solver already initialized")
         return
     end
-    show_info && 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()
@@ -459,7 +459,7 @@ function initialize!(solver::Solver; show_info=true)
         # initialize(problem, ....)
     end
     t1 = round(Base.time()-t0, 2)
-    show_info && info("Initialized solver in $t1 seconds.")
+    info("Initialized solver in $t1 seconds.")
     solver.initialized = true
 end
 
@@ -468,7 +468,7 @@ function get_all_elements(solver::Solver)
     return [elements...;]
 end
 
-function (solver::Solver)(field_name::AbstractString, time::Float64)
+function (solver::Solver)(field_name::String, time::Float64)
     fields = []
     for problem in get_problems(solver)
         field = problem(field_name, time)
@@ -482,14 +482,14 @@ function (solver::Solver)(field_name::AbstractString, time::Float64)
 end
 
 """ Default update for solver. """
-function update!{S}(solver::Solver{S}; show_info=true)
+function update!{S}(solver::Solver{S})
     u = solver.u
     la = solver.la
 
-    show_info && info("Updating problems ...")
+    info("Updating problems ...")
     t0 = Base.time()
 
-    for problem in solver.problems
+    for problem in get_problems(solver)
         assembly = get_assembly(problem)
         elements = get_elements(problem)
         # update solution, first for assembly (u,la) ...
@@ -498,123 +498,45 @@ function update!{S}(solver::Solver{S}; show_info=true)
         update!(problem, assembly, elements, solver.time)
     end
 
-    # if io is attached to solver, update hdf / xml also
-    if !isnull(solver.xdmf)
-        update_xdmf!(solver)
-    end
-
     t1 = round(Base.time()-t0, 2)
-    show_info && info("Updated problems in $t1 seconds.")
+    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)
-
-    # 1. save geometry
-    X_ = solver("geometry", solver.time)
-    node_ids = sort(collect(keys(X_)))
-    X = hcat([X_[nid] for nid in node_ids]...)
-    ndim, nnodes = size(X)
-    geom_type = (ndim == 2 ? "XY" : "XYZ")
-    data_node_ids = new_dataitem(xdmf, "/Node IDs", node_ids)
-    data_geometry = new_dataitem(xdmf, "/Geometry", X)
-    geometry = new_element("Geometry")
-    set_attribute(geometry, "Type", geom_type)
-    add_child(geometry, data_geometry)
-
-    # 2. save topology
-    nid_mapping = Dict(j=>i for (i, j) in enumerate(node_ids))
-    all_elements = get_all_elements(solver)
-    nelements = length(all_elements)
-    debug("Saving topology: $nelements elements total.")
-    element_types = unique(map(get_element_type, all_elements))
-
-    xdmf_element_mapping = Dict(
-        "Poi1" => "Polyvertex",
-        "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")
-
-    topology = []
-    for element_type in element_types
-        elements = filter_by_element_type(element_type, all_elements)
-        nelements = length(elements)
-        info("Xdmf save: $nelements elements of type $element_type")
-        sort!(elements, by=get_element_id)
-        element_ids = map(get_element_id, elements)
-        element_conn = map(element -> [nid_mapping[j]-1 for j in get_connectivity(element)], elements)
-        element_conn = hcat(element_conn...)
-        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_element("Topology")
-        set_attribute(topology_, "TopologyType", xdmf_element_mapping[element_code])
-        set_attribute(topology_, "NumberOfElements", length(elements))
-        add_child(topology_, dataitem)
-        push!(topology, topology_)
-    end
-
-    # 3. save solved field
-    frame = new_element("Grid")
-    time = new_child(frame, "Time")
-    set_attribute(time, "Value", solver.time)
-    add_child(frame, geometry)
-    for topo in topology
-        add_child(frame, topo)
-    end
-
-    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]
+""" Default postprocess for solver. Loop all problems and run postprocess
+functions to calculate secondary fields, i.e. contact pressure, stress,
+heat flux, reaction force etc. quantities.
+"""
+function postprocess!(solver::Solver)
+    info("Running postprocess scripts for solver...")
+    for problem in get_problems(solver)
+        for field_name in problem.postprocess_fields
+            field = Val{Symbol(field_name)}
+            info("Running postprocess for problem $(problem.name), field $field_name")
+            postprocess!(problem, solver.time, field)
         end
-        ndim = 3
     end
-    U = zeros(X)
-    for nid in node_ids
-        loc = nid_mapping[nid]
-        U[:,loc] = U_[nid]
-    end
-    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
-    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, new_dataitem(xdmf, path, U))
-    add_child(frame, attribute)
-    if (S == Linear) || ((S == Nonlinear) && has_converged(solver))
-        add_child(temporal_collection, frame)
-    end
-    save!(xdmf)
 end
 
+""" Default xdmf update for solver. Loop all problems and write them individually
+to Xdmf file. By default write the main unknown field (displacement, temperature,
+...) and any fields requested separately in `problem.postprocess_fields` vector
+(stress, strain, ...)
+"""
+function update_xdmf!(solver::Solver)
+    if isnull(solver.xdmf)
+        info("update_xdmf: xdmf not attached to solver, not writing output to file.")
+        info("turn Xdmf writing on to solver by typing: solver.xdmf = Xdmf(\"results\")")
+        return
+    end
+    xdmf = get(solver.xdmf)
+    for problem in get_problems(solver)
+        fields = [get_unknown_field_name(problem); problem.postprocess_fields]
+        if is_boundary_problem(problem)
+            fields = [fields; get_parent_field_name(problem)]
+        end
+        update_xdmf!(xdmf, problem, solver.time, fields)
+    end
+end
 
 ### Nonlinear quasistatic solver
 
@@ -671,18 +593,23 @@ function (solver::Solver{Nonlinear})()
         info("Increment time t=$(round(solver.time, 3))")
         info(repeat("-", 80))
 
-        # 2.1 update linearized assemblies
+        # 2.1 update assemblies
         assemble!(solver)
+
         # 2.2 call solver for linearized system
         solve!(solver)
+
         # 2.3 update solution back to elements
         update!(solver)
 
         # 2.4 check convergence
-        if has_converged(solver)
+        if properties.iteration >= properties.min_iterations && has_converged(solver)
             info("Converged in $(properties.iteration) iterations.")
-            properties.iteration >= properties.min_iterations && return true
-            info("Convergence criteria met, but iteration < min_iterations, continuing...")
+            # 2.4.1 run any postprocessing of problems
+            postprocess!(solver)
+            # 2.4.2 update Xdmf output
+            update_xdmf!(solver)
+            return true
         end
     end
 
@@ -721,8 +648,8 @@ Main differences in this solver, compared to nonlinear solver are:
 type Linear <: AbstractSolver
 end
 
-function assemble!(solver::Solver{Linear}; show_info=true)
-    show_info && info("Assembling problems ...")
+function assemble!(solver::Solver{Linear})
+    info("Assembling problems ...")
     tic()
     nproblems = 0
     ndofs = 0
@@ -731,27 +658,27 @@ function assemble!(solver::Solver{Linear}; show_info=true)
             assemble!(problem, solver.time)
             nproblems += 1
         else
-            show_info && info("$(problem.name) already assembled, skipping.")
+            info("$(problem.name) already assembled, skipping.")
         end
         ndofs = max(ndofs, size(problem.assembly.K, 2))
     end
     solver.ndofs = ndofs
     t1 = round(toq(), 2)
-    show_info && info("Assembled $nproblems problems in $t1 seconds. ndofs = $ndofs.")
+    info("Assembled $nproblems problems in $t1 seconds. ndofs = $ndofs.")
 end
 
-function (solver::Solver{Linear})(; show_info=true)
+function (solver::Solver{Linear})()
     t0 = Base.time()
-    show_info && info(repeat("-", 80))
-    show_info && info("Starting linear solver")
-    show_info && info("Increment time t=$(round(solver.time, 3))")
-    show_info && info(repeat("-", 80))
+    info(repeat("-", 80))
+    info("Starting linear solver")
+    info("Increment time t=$(round(solver.time, 3))")
+    info(repeat("-", 80))
     initialize!(solver)
     assemble!(solver)
     solve!(solver)
     update!(solver)
     t1 = round(Base.time()-t0, 2)
-    show_info && info("Linear solver ready in $t1 seconds.")
+    info("Linear solver ready in $t1 seconds.")
 end
 
 """ Convenience function to call linear solver. """
@@ -769,60 +696,3 @@ function LinearSolver(name::AbstractString, problems::Problem...)
 end
 
 ### End of linear quasistatic solver
-
-### Postprocessor
-
-type Postprocessor <: AbstractSolver
-    assembly :: Assembly
-    F :: Union{Factorization, Void}
-end
-
-function Postprocessor()
-    Postprocessor(Assembly(), nothing)
-end
-
-function assemble!(solver::Solver{Postprocessor}; show_info=true)
-    show_info && info("Assembling problems ...")
-    tic()
-    nproblems = 0
-    ndofs = 0
-    assembly = solver.properties.assembly
-    empty!(assembly)
-    for problem in get_problems(solver)
-        for element in get_elements(problem)
-            postprocess!(assembly, problem, element, solver.time)
-        end
-        nproblems += 1
-        ndofs = max(ndofs, size(problem.assembly.K, 2))
-    end
-    solver.ndofs = ndofs
-    t1 = round(toq(), 2)
-    show_info && info("Assembled $nproblems problems in $t1 seconds. ndofs = $ndofs.")
-end
-
-function (solver::Solver{Postprocessor})(; show_info=true)
-    t0 = Base.time()
-    show_info && info(repeat("-", 80))
-    show_info && info("Starting postprocessor")
-    show_info && info("Increment time t=$(round(solver.time, 3))")
-    show_info && info(repeat("-", 80))
-    initialize!(solver)
-    assemble!(solver)
-    assembly = solver.properties.assembly
-    M = sparse(assembly.M)
-    f = sparse(assembly.f)
-    F = cholfact(M)
-    q = F \ f
-    t1 = round(Base.time()-t0, 2)
-    show_info && info("Postprocess of results ready in $t1 seconds.")
-    return q
-end
-
-""" Convenience function to call postprocessor. """
-function Postprocessor(problems::Problem...)
-    solver = Solver(Postprocessor, "default postprocessor")
-    if length(problems) != 0
-        push!(solver, problems...)
-    end
-    return solver
-end

src/solvers_modal.jl

@@ -209,8 +209,7 @@ Parameters
 sigma
     Shift stiffness matrix by adding diagonal term, i.e. K_shifted = K + sigma*I
 """
-function (solver::Solver{Modal})(; show_info=true, debug=false,
-              bc_invertible=false, P=nothing, symmetric=true,
+function (solver::Solver{Modal})(; bc_invertible=false, P=nothing, symmetric=true,
               empty_assemblies_before_solution=true, dense=false,
               info_matrices=false, sigma=0.0)
     info(repeat("-", 80))
@@ -273,13 +272,6 @@ function (solver::Solver{Modal})(; show_info=true, debug=false,
 
     info("Calculate $(props.nev) eigenvalues...")
 
-    if debug && length(nz) < 100
-        info("Stiffness matrix:")
-        dump(round(full(K_red)))
-        info("Mass matrix:")
-        dump(round(full(M_red)))
-    end
-
     tic()
  
     if symmetric
@@ -358,14 +350,18 @@ function (solver::Solver{Modal})(; show_info=true, debug=false,
         end
     end
 
-    if !isnull(solver.xdmf)
-        update_xdmf!(solver)
-    end
+    update_xdmf!(solver)
+    
     return true
 
 end
 
-function update_xdmf!(solver::Solver{Modal}; show_info=true)
+function update_xdmf!(solver::Solver{Modal})
+
+    if isnull(solver.xdmf)
+        info("update_xdmf: xdmf not attached to solver, not writing file output.")
+        return
+    end
 
     if maximum(abs(imag(solver.properties.eigvals))) > 1.0e-9
         error("Writing imaginary eigenvalues for Xdmf not supported.")