diff --git a/src/JuliaFEM.jl b/src/JuliaFEM.jl index 95ab2b6..a5b5ff4 100644 --- a/src/JuliaFEM.jl +++ b/src/JuliaFEM.jl @@ -30,7 +30,7 @@ module Testing end include("io.jl") -export Xdmf, h5file, xmffile, has_child, get_child, new_dataitem +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 diff --git a/src/io.jl b/src/io.jl index 0358e68..5ba64fe 100644 --- a/src/io.jl +++ b/src/io.jl @@ -3,100 +3,6 @@ using HDF5 using LightXML -importall LightXML -importall Base - -function new_element(name::String, attrs::Dict) - x = new_element(name) - for (k, v) in attrs - x[k] = v - end - return x -end - -function setindex!(x::XMLElement, content::Any, attr_name::String) - set_attribute(x, attr_name, content) -end - -function haskey(x::XMLElement, key::String) - return has_child(x, key) || has_attribute(x, key) -end - -function has_child(x::XMLElement, child_name::String) - return get_child(x, child_name) != nothing -end - -function get_attribute(x::XMLElement, attr_name::String) - 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::String, attrs::Dict) - 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::String) - '/' in child_name && return nothing - m = match(r"(\w+)\[(.+)\]", child_name) - if m != nothing - child_name = m[1] - end - childs = [] - for child in child_elements(x) - if name(child) == child_name - push!(childs, child) - end - end - length(childs) == 0 && return nothing - m == nothing && return first(childs) - j = tryparse(Int, m[2]) - isnull(j) || return childs[get(j)] - m[2] == "end" && return childs[end] - m2 = match(r"@(.+)=(.+)", m[2]) - m2 == nothing && throw("Unable to parse: $(m[2])") - attr_name = convert(String, m2[1]) - attr_value = convert(String, m2[2]) - for child in childs - has_attribute(child, attr_name) || continue - if get_attribute(child, attr_name) == attr_value - return child - end - end - throw("Unable to parse: $(m[2])") -end - -function getindex(x::XMLElement, attr_name::String) - attr_name = strip(attr_name, '/') - child = get_child(x, attr_name) - child == nothing || return child - has_attribute(x, attr_name) && return get_attribute(x, attr_name) - if '/' in attr_name - items = map(String, split(attr_name, '/')) - attr_name = first(items) - length(items) > 1 || throw(KeyError(attr_name)) - haskey(x, attr_name) || throw(KeyError(attr_name)) - path = join(items[2:end], '/') - new_item = getindex(x, attr_name) - return new_item[path] - else - throw(KeyError(attr_name)) - end -end type Xdmf name :: String @@ -104,24 +10,6 @@ type Xdmf hdf :: HDF5File end -function new_child(xdmf::Xdmf, args...; kwargs...) - new_child(xdmf.xml, args...; kwargs...) -end - -function read(xdmf::Xdmf, path::String) - result = getindex(xdmf.xml, path) - if endswith(path, "DataItem") - format = get_attribute(result, "Format") - @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) - else - return result - end -end - function Xdmf() return Xdmf(tempname()) end @@ -134,6 +22,7 @@ function xmffile(xdmf::Xdmf) return xdmf.name*".xmf" end +""" Initialize a new Xdmf object. """ function Xdmf(name::String; overwrite=false) xdmf = new_element("Xdmf") h5file = "$name.h5" @@ -164,6 +53,175 @@ function Xdmf(name::String; overwrite=false) return Xdmf(name, xdmf, 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 + + + + + + + + + +""" +function get_temporal_collection(xdmf::Xdmf) + domain = find_element(xdmf.xml, "Domain") + grid = nothing + if domain == nothing + debug("Xdmf: creating new temporal collection") + domain = new_child(xdmf.xml, "Domain") + grid = new_child(domain, "Grid") + set_attribute(grid, "CollectionType", "Temporal") + set_attribute(grid, "GridType", "Collection") + end + grid = find_element(domain, "Grid") + return grid +end + +""" Returns some spesific child xml element from a array of XMLElement based on, +"Xdmf extensions" see [1] for details. + +Parameters +---------- +child_elements :: Vector{XMLElement} + A vector of XMLElements where to perform filtering. +child_name :: String + Child element name, maybe containing Xdmf instructions + +Returns +------- +nothing if nothing is found, otherwise XMLElement matching to filtering + +Examples +-------- +julia> grid1 = new_element("Grid") +julia> add_text(grid1, "I am first grid") +julia> grid2 = new_element("Grid") +julia> add_text(grid2, "I am second grid") +julia> set_attribute(grid2, "Name", "Frame 2") +julia> grid3 = new_element("Grid") +julia> add_text(grid3, "I am third grid") +julia> grids = [grid1, grid2, grid3] + +To return second Grid element, one can use + +julia> xdmf_filter(grids, "Grid[2]") + +To return Grid which has attribute Name="Frame 2", use + +julia> xdmf_filter(grids, "Grid[@name=Frame 2]") + +To pick last Grid, use [end], e.g. + +julia> xdmf_filter(grids, "Grid[end]"). + +References +---------- +[1] http://www.xdmf.org/index.php/XDMF_Model_and_Format +""" +function xdmf_filter(child_elements, child_name) + if '/' in child_name # needs path traversal + return nothing + end + + # filter children elements using syntax child[X] -> rename child_name + m = match(r"(\w+)\[(.+)\]", child_name) + if m != nothing + child_name = m[1] + end + + # first find any relevant child elements (has same tag) + childs = [] + for child in child_elements + if LightXML.name(child) == child_name + push!(childs, child) + end + end + + # childs not found at all + length(childs) == 0 && return nothing + + # by default return first + m == nothing && return first(childs) + + # if [end] return last + m[2] == "end" && return childs[end] + + # otherwise try parse int and return nth children from list + parsed_int = tryparse(Int, m[2]) + if !isnull(parsed_int) + idx = get(parsed_int) + if (idx > 0) && (idx <= length(childs)) + return childs[idx] + else + # wrong index + return nothing + end + end + + # [X] is something else than integer, filter children elements using syntax child[@attr=value] + m2 = match(r"@(.+)=(.+)", m[2]) + m2 == nothing && throw("Unable to parse: $(m[2])") + attr_name = convert(String, m2[1]) + attr_value = convert(String, m2[2]) + for child in childs + has_attribute(child, attr_name) || continue + if attribute(child, attr_name) == attr_value + return child + end + end + + # nothing found + return nothing +end + +""" Traverse XML path. Xdmf filtering can be used, so it's possible to find +data from xml using syntax e.g. + +julia> traverse(xdmf, x, "/Domain/Grid[2]/Grid[@Name=Frame 1]/DataItem") +""" +function traverse(xdmf::Xdmf, x::XMLElement, attr_name::String) + attr_name = strip(attr_name, '/') + + if has_attribute(x, attr_name) + return attribute(x, attr_name) + end + + childs = child_elements(x) + + if '/' in attr_name + items = split(attr_name, '/') + new_item = xdmf_filter(childs, first(items)) + new_path = join(items[2:end], '/') + return traverse(xdmf, new_item, new_path) + end + + child = xdmf_filter(childs, attr_name) + return child +end + +""" Read data from Xdmf file. + +Traversing is supported, so one can easily traverse XML tree e.g. +julia> read(xdmf, "/Domain/Grid/Grid[2]/Geometry") +""" +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) + else + return result + end +end + + function save!(xdmf::Xdmf) doc = XMLDocument() set_root(doc, xdmf.xml) diff --git a/src/solvers.jl b/src/solvers.jl index 095e7c7..79b5709 100644 --- a/src/solvers.jl +++ b/src/solvers.jl @@ -488,20 +488,6 @@ function (solver::Solver)(field_name::AbstractString, time::Float64) return merge(fields...) end -function get_temporal_collection(xdmf::Xdmf) - domain = find_element(xdmf.xml, "Domain") - grid = nothing - if domain == nothing - info("Xdmf: creating new temporal collection") - domain = new_child(xdmf.xml, "Domain") - grid = new_child(domain, "Grid") - set_attribute(grid, "CollectionType", "Temporal") - set_attribute(grid, "GridType", "Collection") - end - grid = find_element(domain, "Grid") - return grid -end - """ Default update for solver. """ function update!{S}(solver::Solver{S}; show_info=true) u = solver.u @@ -540,7 +526,8 @@ function update_xdmf!{S}(solver::Solver{S}; show_info=true) 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", Dict("Type" => geom_type)) + geometry = new_element("Geometry") + set_attribute(geometry, "Type", geom_type) add_child(geometry, data_geometry) # 2. save topology @@ -588,7 +575,8 @@ function update_xdmf!{S}(solver::Solver{S}; show_info=true) # 3. save solved field frame = new_element("Grid") - new_child(frame, "Time", Dict("Value" => solver.time)) + time = new_child(frame, "Time") + set_attribute(time, "Value", solver.time) add_child(frame, geometry) for topo in topology add_child(frame, topo) diff --git a/src/solvers_modal.jl b/src/solvers_modal.jl index 5a4c870..db7923b 100644 --- a/src/solvers_modal.jl +++ b/src/solvers_modal.jl @@ -420,10 +420,12 @@ function update_xdmf!(solver::Solver{Modal}; show_info=true) info("Creating frequency frame f=$(round(freq, 3)), path=$path") frame = new_element("Grid") - new_child(frame, "Time", Dict("Value" => freq)) + time = new_child(frame, "Time") + set_attribute(time, "Value", freq) # add geometry - geometry = new_element("Geometry", Dict("Type" => geom_type)) + geometry = new_element("Geometry") + set_attribute(geometry, "Type", geom_type) data_node_ids = new_dataitem(xdmf, "/Node IDs", node_ids) data_geometry = new_dataitem(xdmf, "/Geometry", X) add_child(geometry, data_geometry) diff --git a/test/test_io.jl b/test/test_io.jl index 69da349..13537fb 100644 --- a/test/test_io.jl +++ b/test/test_io.jl @@ -23,25 +23,48 @@ end @test attribute(obj2, "Name") == "Test Domain" end -@testset "put and get to xdmf" begin +@testset "Xdmf filtering" begin + grid1 = new_element("Grid") + add_text(grid1, "I am first grid") + grid2 = new_element("Grid") + add_text(grid2, "I am second grid") + set_attribute(grid2, "Name", "Frame 2") + grid3 = new_element("Grid") + add_text(grid3, "I am third grid") + grids = [grid1, grid2, grid3] + + @test content(xdmf_filter(grids, "Grid")) == "I am first grid" + @test content(xdmf_filter(grids, "Grid[1]")) == "I am first grid" + @test content(xdmf_filter(grids, "Grid[2]")) == "I am second grid" + @test content(xdmf_filter(grids, "Grid[3]")) == "I am third grid" + @test content(xdmf_filter(grids, "Grid[end]")) == "I am third grid" + @test content(xdmf_filter(grids, "Grid[@Name=Frame 2]")) == "I am second grid" + @test xdmf_filter(grids, "Grid[0]") == nothing + @test xdmf_filter(grids, "Grid[4]") == nothing + @test xdmf_filter(grids, "Grid[@Name=Frame 3]") == nothing + @test xdmf_filter(grids, "Domain/Grid[@Name=Frame 3]") == nothing + @test xdmf_filter(grids, "Domain") == nothing +end + +@testset "XML traverse" begin xdmf = Xdmf() - domain = new_child(xdmf, "Domain") + domain = new_child(xdmf.xml, "Domain") grid = new_child(domain, "Grid") - grid["CollectionType"] = "Temporal" - grid["GridType"] = "Collection" + set_attribute(grid, "CollectionType", "Temporal") + set_attribute(grid, "GridType", "Collection") frame1 = new_child(grid, "Grid") time1 = new_child(frame1, "Time") - set_attributes(time1, Dict("Value" => 0.0)) + set_attribute(time1, "Value", 0.0) X1 = new_child(frame1, "Geometry") - set_attributes(X1, Dict("Type" => "XY")) + set_attribute(X1, "Type", "XY") frame2 = new_child(grid, "Grid") - set_attributes(frame2, Dict("Name" => "Frame 2")) + set_attribute(frame2, "Name", "Frame 2") time2 = new_child(frame2, "Time") - set_attributes(time2, Dict("Value" => 1.0)) + set_attribute(time2, "Value", 1.0) X2 = new_child(frame2, "Geometry") - set_attributes(X2, Dict("Type" => "XY")) + set_attribute(X2, "Type", "XY") add_child(grid, frame1) add_child(grid, frame2) @@ -50,178 +73,9 @@ end add_child(X2, dataitem) println(xdmf.xml) - @test has_child(xdmf.xml, "Domain") - @test isa(get_child(xdmf.xml, "Domain"), XMLElement) - @test !has_attribute(xdmf.xml, "Domain") - @test isapprox(read(xdmf, "/Domain/Grid/Grid/Time/Value"), 0.0) - @test isapprox(read(xdmf, "/Domain/Grid/Grid[2]/Time/Value"), 1.0) - @test isapprox(read(xdmf, "/Domain/Grid/Grid[end]/Time/Value"), 1.0) - @test isapprox(read(xdmf, "/Domain/Grid/Grid[@Name=Frame 2]/Time/Value"), 1.0) + @test read(xdmf, "/Domain/Grid/Grid/Time/Value") == "0.0" + @test read(xdmf, "/Domain/Grid/Grid[2]/Time/Value") == "1.0" + @test read(xdmf, "/Domain/Grid/Grid[end]/Time/Value") == "1.0" + @test read(xdmf, "/Domain/Grid/Grid[@Name=Frame 2]/Time/Value") == "1.0" @test isapprox(read(xdmf, "/Domain/Grid/Grid[2]/Geometry/DataItem"), [1.0, 2.0]) end - - -#= - -@testset "higher level xdmf" begin - e1 = Element(Quad4, 1, [1, 2, 3, 4]) - e2 = Element(Quad4, 2, [5, 6, 7, 8]) - p1 = Problem(Elasticity, "Body 1", 2) - p2 = Problem(Elasticity, "Body 2", 2) - push!(p1, e1) - push!(p2, e2) - #update!(p1) - - e3 = Element(Seg2, 3, [1, 2]) - e4 = Element(Seg2, 4, [3, 4]) - e5 = Element(Seg2, 5, [5, 6]) - p3 = Problem(Dirichlet, "Fixed BC", 2, "displacement") - p4 = Problem(Contact, "Contact between bodies 1 and 2", 2, "displacement") - push!(p3, e3) - push!(p4, e4) - X = Dict{Int64, Vector{Float64}}( - 1 => [0.0, 0.0], - 2 => [1.0, 0.0], - 3 => [1.0, 1.0], - 4 => [0.0, 1.0], - 5 => [0.0, 2.0], - 6 => [1.0, 2.0], - 7 => [1.0, 3.0], - 8 => [0.0, 3.0]) - u = Dict{Int64, Vector{Float64}}( - 1 => [0.1, 0.1], - 2 => [0.1, 0.1], - 3 => [0.1, 0.1], - 4 => [0.1, 0.1], - 5 => [0.1, 0.1], - 6 => [0.1, 0.1], - 7 => [0.1, 0.1], - 8 => [0.1, 0.1]) - n = Dict{Int64, Vector{Float64}}( - 3 => [0.0, 1.0], - 4 => [0.0, 1.0]) - R = Dict{Int64, Vector{Float64}}( - 1 => [0.0, 1.0], - 2 => [0.0, 1.0]) - update!(e4, "master elements", [e3]) -end - -@testset "save results to disk, linear solver" 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", 0.0 => 12.0) - update!(element, "temperature load", 1.0 => 24.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) - xdmf = Xdmf() - solver = Solver(Linear, problem, bc) - solver.xdmf = xdmf - - solver.time = 0.0 - solver() - empty!(problem.assembly) - solver.time = 1.0 - solver() - - info(solver("temperature", 0.0)) - info(solver("temperature", 1.0)) - info(element("temperature load", [0.0, 0.0], 0.0)) - info(element("temperature load", [0.0, 0.0], 1.0)) - - info("h5 file = $(h5file(xdmf))") - E = read(xdmf.hdf, "/Topology/Quad4/Element IDs") - C = read(xdmf.hdf, "/Topology/Quad4/Connectivity") - N = read(xdmf.hdf, "/Node IDs") - X = read(xdmf.hdf, "/Geometry") - T1 = read(xdmf.hdf, "/Results/Time 0.0/Nodal Fields/Temperature") - T2 = read(xdmf.hdf, "/Results/Time 1.0/Nodal Fields/Temperature") - @test isapprox(E, [-1]) - @test isapprox(C, [0 1 2 3]) - @test isapprox(N, [1, 2, 3, 4]) - X_expected = [0.0 0.0; 1.0 0.0; 1.0 1.0; 0.0 1.0]' - T1_expected = [0.0 0.0 1.0 1.0] - T2_expected = [0.0 0.0 2.0 2.0] - @test isapprox(X, X_expected) - @test isapprox(T1, T1_expected) - @test isapprox(T2, T2_expected) - @test isapprox(read(xdmf, "/Domain/Grid/Grid/Time/Value"), 0.0) - @test read(xdmf, "/Domain/Grid/Grid/Geometry/Type") == "XY" - @test isapprox(read(xdmf, "/Domain/Grid/Grid/Geometry/DataItem"), X_expected) - @test isapprox(read(xdmf, "/Domain/Grid/Grid/Topology/DataItem"), [0 1 2 3]) - @test isapprox(read(xdmf, "/Domain/Grid/Grid/Topology[@TopologyType=Polyline]/DataItem"), [0 1]) - @test isapprox(read(xdmf, "/Domain/Grid/Grid[1]/Attribute[@Name=Temperature]/DataItem"), T1_expected) - @test isapprox(read(xdmf, "/Domain/Grid/Grid[2]/Attribute[@Name=Temperature]/DataItem"), T2_expected) - @test isapprox(read(xdmf, "/Domain/Grid/Grid[end]/Time/Value"), 1.0) - @test isapprox(read(xdmf, "/Domain/Grid/Grid[end]/Topology/DataItem"), [0 1 2 3]) -end - -@testset "save results to disk, nonlinear solver" 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", 0.0 => 12.0) - update!(element, "temperature load", 1.0 => 24.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(Nonlinear, problem, bc) - solver.xdmf = Xdmf() - - solver.time = 0.0 - solver() - solver.time = 1.0 - solver() - - xdmf = get(solver.xdmf) - info("h5 file = $(h5file(xdmf))") - E = read(xdmf.hdf, "/Topology/Quad4/Element IDs") - C = read(xdmf.hdf, "/Topology/Quad4/Connectivity") - N = read(xdmf.hdf, "/Node IDs") - X = read(xdmf.hdf, "/Geometry") - T11 = read(xdmf.hdf, "/Results/Time 0.0/Iteration 1/Nodal Fields/Temperature") - T12 = read(xdmf.hdf, "/Results/Time 0.0/Iteration 2/Nodal Fields/Temperature") - T21 = read(xdmf.hdf, "/Results/Time 1.0/Iteration 1/Nodal Fields/Temperature") - T22 = read(xdmf.hdf, "/Results/Time 1.0/Iteration 2/Nodal Fields/Temperature") - X_expected = [ - 0.0 0.0 - 1.0 0.0 - 1.0 1.0 - 0.0 1.0] - T1_expected = [0.0 0.0 1.0 1.0] - T2_expected = [0.0 0.0 2.0 2.0] - @test isapprox(T12, T1_expected) - @test isapprox(T22, T2_expected) - @test isapprox(read(xdmf, "/Domain/Grid/Grid/Time/Value"), 0.0) - @test read(xdmf, "/Domain/Grid/Grid/Geometry/Type") == "XY" - @test isapprox(read(xdmf, "/Domain/Grid/Grid/Geometry/DataItem"), X_expected') - @test isapprox(read(xdmf, "/Domain/Grid/Grid/Topology/DataItem"), [0 1 2 3]) - @test isapprox(read(xdmf, "/Domain/Grid/Grid/Topology[@TopologyType=Polyline]/DataItem"), [0 1]) - @test isapprox(read(xdmf, "/Domain/Grid/Grid[1]/Attribute[@Name=Temperature]/DataItem"), T1_expected) - @test isapprox(read(xdmf, "/Domain/Grid/Grid[2]/Attribute[@Name=Temperature]/DataItem"), T2_expected) - @test isapprox(read(xdmf, "/Domain/Grid/Grid[end]/Time/Value"), 1.0) - @test isapprox(read(xdmf, "/Domain/Grid/Grid[end]/Topology/DataItem"), [0 1 2 3]) -end - -=#