io

src/JuliaFEM.jl

@@ -8,6 +8,21 @@ module JuliaFEM
 
 importall Base
 
+module Testing
+
+    if VERSION >= v"0.5-"
+        using Base.Test
+    else
+        using BaseTestNext
+    end
+
+    export @test, @testset, @test_throws
+
+end
+
+include("io.jl")
+export ModelIO
+
 include("fields.jl")
 export Field, DCTI, DVTI, DCTV, DVTV, CCTI, CVTI, CCTV, CVTV, Increment
 include("types.jl")  # data types: Point, IntegrationPoint, ...
@@ -15,7 +30,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
+export Node, AbstractElement, Element, update!, get_connectivity, get_basis,
+       get_dbasis, inside, get_local_coordinates
 include("elements_lagrange.jl") # Continuous Galerkin (Lagrange) elements
 export get_reference_coordinates, get_interpolation_polynomial
 export Poi1,
@@ -35,7 +51,7 @@ include("integrate.jl")  # default integration points for elements
 export get_integration_points
 
 include("sparse.jl")
-export add!, SparseMatrixCOO, get_nonzero_rows
+export add!, SparseMatrixCOO, get_nonzero_rows, get_nonzero_columns
 
 include("problems.jl") # common problem routines
 export Problem, AbstractProblem, FieldProblem, BoundaryProblem,
@@ -78,11 +94,8 @@ include("problems_mortar.jl")
 include("problems_mortar_2d.jl")
 include("problems_mortar_3d.jl")
 include("problems_mortar_2d_autodiff.jl")
-export calculate_normals,
-       calculate_normals!,
-       project_from_slave_to_master,
-       project_from_master_to_slave,
-       Mortar, get_slave_elements,
+export calculate_normals, calculate_normals!, project_from_slave_to_master,
+       project_from_master_to_slave, Mortar, get_slave_elements,
        get_polygon_clip
 
 ### Mortar methods, contact mechanics extension ###
@@ -92,22 +105,11 @@ include("problems_contact_3d.jl")
 include("problems_contact_2d_autodiff.jl")
 export Contact
 
-#=
-module API
-include("api.jl")
-# export ....
-end
-=#
-
 module Preprocess
 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!,
-       find_nearest_nodes,
-       reorder_element_connectivity!
+export create_elements, Mesh, add_node!, add_nodes!, add_element!,
+       add_elements!, add_element_to_element_set!, add_node_to_node_set!,
+       find_nearest_nodes, reorder_element_connectivity!
 include("preprocess_abaqus_reader.jl")
 export parse_abaqus, parse_section, parse_element_section,
        abaqus_read_mesh, abaqus_read_model
@@ -130,50 +132,21 @@ end
 export get_mesh, get_model
 
 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
+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!
-export DataFrame
+export XDMF, xdmf_new_result!, xdmf_save_field!, xdmf_save!, DataFrame
+
 end
 export Postprocessor
 
-# This connects model from preprocess_abaqus_reader to 
-# other JuliaFEM ecosystem and solves problem.
 module Abaqus
 include("abaqus.jl")
 export abaqus_read_model, abaqus_run_model, abaqus_open_results
 end
 
-""" JuliaFEM testing routines. """
-module Testing
-if VERSION >= v"0.5-"
-    using Base.Test
-else
-    using BaseTestNext
-end
-
-export @test, @testset, @test_throws
-#include("test.jl")
-end
-
-#=
-module MaterialModels
-include("vonmises.jl")
-end
-=#
-
-#=
-module Interfaces
-include("interfaces.jl")
-end
-=#
-
-end # module
+end #

src/io.jl

@@ -0,0 +1,53 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
+
+using HDF5
+using LightXML
+
+type ModelIO
+    name :: AbstractString
+    xdmf :: XMLElement
+end
+
+function ModelIO()
+    return ModelIO(tempname())
+end
+
+function ModelIO(name::AbstractString)
+    xdmf = new_element("Xdmf")
+    set_attribute(xdmf, "xmlns:xi", "http://www.w3.org/2001/XInclude")
+    set_attribute(xdmf, "Version", "2.1")
+    return ModelIO(name, xdmf)
+end
+
+function h5file(mio::ModelIO)
+    return mio.name*".h5"
+end
+
+function put!{T,N}(mio::ModelIO, path::AbstractString, data::Array{T,N})
+    hdf = h5file(mio)
+    h5write(hdf, path, data)
+end
+
+function get_dataitem{T,N}(mio::ModelIO, path::AbstractString, data::Array{T,N}; format="HDF")
+    dataitem = new_element("DataItem")
+    n, m = size(data)
+    set_attribute(dataitem, "DataType", "$T")
+    set_attribute(dataitem, "Dimensions", "$n $m")
+    set_attribute(dataitem, "Format", format)
+    if format == "HDF"
+        hdf = basename(h5file(mio))
+        add_text(dataitem, "$hdf:$path")
+    end
+    return dataitem
+end
+
+function get(mio::ModelIO, path::AbstractString)
+    h5read(mio.name*".h5", path)
+end
+
+function save!(mio::ModelIO)
+    doc = XMLDocument()
+    set_root(doc, mio.xdmf)
+    save_file(doc, mio.name*".xmf")
+end

src/postprocess_utils.jl

@@ -139,7 +139,7 @@ function calc_nodal_values!(problem::Problem, field_name::AbstractString, field_
 end
 
 """
-Return node ids + vector of values 
+Return node ids + vector of values
 """
 function get_nodal_vector(elements::Vector, field_name::AbstractString, time::Float64)
     f = Dict()
@@ -201,17 +201,22 @@ end
 
 """ Return field calculated to nodal points for elements in problem p. """
 function call(problem::Problem, field_name::AbstractString, time::Float64=0.0)
-    f = Dict()
+    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
-            f[c] = v
         end
     end
     return f
@@ -428,4 +433,3 @@ function getindex(problem::Problem, field_name::AbstractString)
     end
     return DVTV(increments)
 end
-

src/solvers.jl

@@ -4,24 +4,25 @@
 abstract AbstractSolver
 
 type Solver{S<:AbstractSolver}
-    name :: AbstractString # some descriptive name for problem
-    time :: Float64        # current time
+    name :: AbstractString       # some descriptive name for problem
+    time :: Float64              # current time
     problems :: Vector{Problem}
-    norms :: Vector{Tuple}  # solution norms for convergence studies
-    ndofs :: Int        # number of degrees of freedom in problem
+    norms :: Vector{Tuple}       # solution norms for convergence studies
+    ndofs :: Int                 # number of degrees of freedom in problem
+    io :: Nullable{ModelIO}      # input/output handle
     properties :: S
 end
 
 
 function Solver{S<:AbstractSolver}(::Type{S}, name="solver", properties...)
     variant = S(properties...)
-    solver = Solver{S}(name, 0.0, [], [], 0, variant)
+    solver = Solver{S}(name, 0.0, [], [], 0, nothing, variant)
     return solver
 end
 
 function Solver{S<:AbstractSolver}(::Type{S}, problems::Problem...)
     variant = S()
-    solver = Solver{S}("$(S)Solver", 0.0, [], [], 0, variant)
+    solver = Solver{S}("$(S)Solver", 0.0, [], [], 0, nothing, variant)
     push!(solver.problems, problems...)
     return solver
 end
@@ -430,6 +431,44 @@ function initialize!(solver::Solver; show_info=true)
     show_info && info("Initialized problems in $t1 seconds.")
 end
 
+function get_all_elements(solver::Solver)
+    elements = [get_elements(problem) for problem in get_problems(solver)]
+    return [elements...;]
+end
+
+function get_element_type{E}(element::Element{E})
+    return E
+end
+
+function get_element_id{E}(element::Element{E})
+    return element.id
+end
+
+function is_element_type{E}(element::Element{E}, element_type)
+    return is(E, element_type)
+end
+
+function filter_by_element_type(element_type, elements)
+    return filter(element -> is_element_type(element, element_type), elements)
+end
+
+function call(solver::Solver, field_name::AbstractString, time::Float64)
+    fields = [problem(field_name, time) for problem in get_problems(solver)]
+    return merge(fields...)
+end
+
+function get_temporal_collection(mio::ModelIO)
+    grid = find_element(mio.xdmf, "Grid")
+    if grid == nothing
+        info("Xdmf: creating new temporal collection")
+        domain = new_child(mio.xdmf, "Domain")
+        grid = new_child(domain, "Grid")
+        set_attribute(grid, "CollectionType", "Temporal")
+        set_attribute(grid, "GridType", "Collection")
+    end
+    return grid
+end
+
 """ Default update for solver. """
 function update!(solver::Solver, u::Vector, la::Vector; show_info=true)
     show_info && info("Updating problems ...")
@@ -442,6 +481,86 @@ function update!(solver::Solver, u::Vector, la::Vector; show_info=true)
         # .. and then from assembly (u,la) to elements
         update!(problem, assembly, elements, solver.time)
     end
+
+    # if io is attached to solver, update hdf / xml also
+    if !isnull(solver.io)
+        io = get(solver.io)
+        xdmf = io.xdmf
+        temporal_collection = get_temporal_collection(io)
+        frame = new_child(temporal_collection, "Grid")
+        time_item = new_child(frame, "Time")
+        set_attribute(time_item, "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]...)
+        put!(io, "/Node IDs", node_ids)
+        path = "/Geometry"
+        put!(io, path, geometry)
+        dataitem = get_dataitem(io, path, geometry)
+
+        ndim, nnodes = size(geometry)
+        geom_type = ndim == 2 ? "XY" : "XYZ"
+        geom = new_child(frame, "Geometry")
+        set_attribute(geom, "Type", geom_type)
+        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))
+
+        element_mapping = Dict(
+            "Quad4" => "Quadrilateral",
+            "Seg2" => "Polyline",
+            )
+
+        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]
+            put!(io, "/Topology/$element_code/Element IDs", element_ids)
+            path = "/Topology/$element_code/Connectivity"
+            put!(io, path, element_conn)
+            dataitem = get_dataitem(io, path, element_conn)
+            topology = new_child(frame, "Topology")
+            set_attribute(topology, "TopologyType", 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)
+        path = "/Results/$time/Nodal Fields/$unknown_field_name"
+        put!(io, path, U)
+        dataitem = get_dataitem(io, path, U)
+        attribute = new_child(frame, "Attribute")
+        set_attribute(attribute, "Name", unknown_field_name)
+        set_attribute(attribute, "Center", field_center)
+        add_child(attribute, dataitem)
+        save!(io)
+    end
+
     t1 = round(Base.time()-t0, 2)
     show_info && info("Updated problems in $t1 seconds.")
 end
@@ -685,4 +804,3 @@ function Postprocessor(name::AbstractString, problems::Problem...)
     solver.name = name
     return solver
 end
-

test/test_io.jl

@@ -0,0 +1,53 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
+
+using JuliaFEM
+using JuliaFEM.Testing
+importall Base
+
+@testset "create new result" begin
+    r = ModelIO()
+    expected = """<Xdmf xmlns:xi="http://www.w3.org/2001/XInclude" Version="2.1"/>"""
+    @test string(r.xdmf) == expected
+end
+
+@testset "put and get result" begin
+    r = ModelIO()
+    put!(r, "/1/2/3", [1 2 3])
+    @test isapprox(get(r, "/1/2/3"), [1 2 3])
+end
+
+@testset "save results to disk" begin
+    X = Dict{Int, Vector{Float64}}(
+        1 => [0.0,0.0],
+        2 => [1.0,0.0],
+        3 => [1.0,1.0],
+        4 => [0.0,1.0])
+    element = Element(Quad4, [1, 2, 3, 4])
+    update!(element, "geometry", X)
+    update!(element, "temperature thermal conductivity", 6.0)
+    update!(element, "temperature load", 12.0)
+    problem = Problem(Heat, "one element heat problem", 1)
+    problem.properties.formulation = "2D"
+    push!(problem, element)
+    boundary_element = Element(Seg2, [1, 2])
+    update!(boundary_element, "geometry", X)
+    update!(boundary_element, "temperature 1", 0.0)
+    bc = Problem(Dirichlet, "fixed", 1, "temperature")
+    push!(bc, boundary_element)
+    solver = Solver(Linear, problem, bc)
+    solver.io = ModelIO()
+    solver()
+    io = get(solver.io)
+    info("h5 file = $(io.name).h5")
+    E = get(io, "/Topology/Quad4/Element IDs")
+    C = get(io, "/Topology/Quad4/Connectivity")
+    N = get(io, "/Node IDs")
+    X = get(io, "/Geometry")
+    T = get(io, "/Results/Time 0.0/Nodal Fields/Temperature")
+    @test isapprox(E, [-1])
+    @test isapprox(C, [0 1 2 3])
+    @test isapprox(N, [1, 2, 3, 4])
+    @test isapprox(X, [0.0 0.0; 1.0 0.0; 1.0 1.0; 0.0 1.0]')
+    @test isapprox(T, [0.0 0.0 1.0 1.0])
+end