diff --git a/src/JuliaFEM.jl b/src/JuliaFEM.jl index 22b1cba..c8d0f36 100644 --- a/src/JuliaFEM.jl +++ b/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 # diff --git a/src/io.jl b/src/io.jl new file mode 100644 index 0000000..977366f --- /dev/null +++ b/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 diff --git a/src/postprocess_utils.jl b/src/postprocess_utils.jl index f6b169b..70d51b3 100644 --- a/src/postprocess_utils.jl +++ b/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 - diff --git a/src/solvers.jl b/src/solvers.jl index ab8431a..5fd06d7 100644 --- a/src/solvers.jl +++ b/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 - diff --git a/test/test_io.jl b/test/test_io.jl new file mode 100644 index 0000000..b3ce12e --- /dev/null +++ b/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 = """""" + @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