diff --git a/src/JuliaFEM.jl b/src/JuliaFEM.jl index d62e219..b5c1353 100644 --- a/src/JuliaFEM.jl +++ b/src/JuliaFEM.jl @@ -138,6 +138,7 @@ export calc_nodal_values!, calculate_second_moment_of_mass include("postprocess_xdmf.jl") export XDMF, xdmf_new_result!, xdmf_save_field!, xdmf_save! +export DataFrame end export Postprocessor @@ -145,7 +146,7 @@ export Postprocessor # other JuliaFEM ecosystem and solves problem. module Abaqus include("abaqus.jl") -export abaqus_read_model, abaqus_run_model, abaqus_read_results, abaqus_run_test +export abaqus_read_model, abaqus_run_model, abaqus_open_results end """ JuliaFEM testing routines. """ diff --git a/src/abaqus.jl b/src/abaqus.jl index e06a97f..b2da1dc 100644 --- a/src/abaqus.jl +++ b/src/abaqus.jl @@ -6,6 +6,7 @@ importall Base using JuliaFEM using JuliaFEM.Preprocess using JuliaFEM.Postprocess +using LightXML ### Model definitions for ABAQUS data model @@ -17,6 +18,8 @@ abstract AbstractBoundaryCondition abstract AbstractOutputRequest type Model + path :: AbstractString + name :: AbstractString mesh :: Mesh materials :: Dict{Symbol, AbstractMaterial} properties :: Vector{AbstractProperty} @@ -211,7 +214,9 @@ end function abaqus_read_model(fn; read_mesh=true) - model = Model(Mesh(), Dict(), Vector(), Vector(), Vector(), Vector()) + model_path = dirname(fn) + model_name = first(splitext(basename(fn))) + model = Model(model_path, model_name, Mesh(), Dict(), [], [], [], []) if read_mesh model.mesh = abaqus_read_mesh(fn) @@ -230,7 +235,6 @@ function abaqus_read_model(fn; read_mesh=true) end end close(fid) - maybe_close_section!(model, state) return model @@ -503,27 +507,27 @@ function get_child_element(element_type::Symbol, element_side::Symbol, :P6 => (:Quad4, [4, 8, 5, 1])) ) - if !haskey(element_mapping, element_type) - error("Unable to find child element for element of type $element_type for side $element_side, check mapping.") + if !haskey(element_mapping, element_type) + error("Unable to find child element for element of type $element_type for side $element_side, check mapping.") + end + + if !haskey(element_mapping[element_type], element_side) + error("Unable to find child element side mapping for element of type $element_type for side $element_side, check mapping.") + end + + child_element, child_element_lconn = element_mapping[element_type][element_side] + child_element_gconn = element_connectivity[child_element_lconn] + return child_element, child_element_lconn, child_element_gconn end - if !haskey(element_mapping[element_type], element_side) - error("Unable to find child element side mapping for element of type $element_type for side $element_side, check mapping.") + function determine_solver_type(model::Model, step::AbstractStep) + # FIXME + return Linear end - child_element, child_element_lconn = element_mapping[element_type][element_side] - child_element_gconn = element_connectivity[child_element_lconn] - return child_element, child_element_lconn, child_element_gconn -end - -function determine_solver_type(model::Model, step::AbstractStep) - # FIXME - return Linear -end - -function process_output_request(model::Model, solver::Solver, output_request::AbstractOutputRequest) - kind = Val{output_request.kind} - target = Val{output_request.target} + function process_output_request(model::Model, solver::Solver, output_request::AbstractOutputRequest) + kind = Val{output_request.kind} + target = Val{output_request.target} process_output_request(model, solver, output_request, kind, target) end @@ -534,20 +538,22 @@ function process_output_request(model::Model, solver::Solver, output_request::Ab data = output_request.data options = output_request.options info("nodal output request with data $data and options $options") + code_mapping = Dict(:U => "displacement", :COORD => "geometry") + abbr_mapping = Dict(:U => :U, :COORD => :COOR) # ..? for row in data + tables = Any[] for code in row + haskey(code_mapping, code) || continue for problem in model.problems - if code == :U - vals = problem("displacement", solver.time) - node_ids = sort(collect(keys(vals))) - u1 = [vals[id][1] for id in node_ids] - u2 = [vals[id][2] for id in node_ids] - u3 = [vals[id][3] for id in node_ids] - d = DataFrame(; id=node_ids, u1=u1, u2=u2, u3=u3) - println(d) - end + field_name = code_mapping[code] + abbr = abbr_mapping[code] + table = problem(DataFrame, field_name, abbr, solver.time) + push!(tables, table) end end + length(tables) != 0 || continue + results = join(tables..., on=:id, kind=:outer) + println(results) end end @@ -590,53 +596,76 @@ function call(model::Model) end end - return true + return 0 end -function abaqus_read_results +function abaqus_download(name) + fn = "$name.inp" + if !haskey(ENV, "ABAQUS_DOWNLOAD_URL") + info(""" + ABAQUS input file $fn not found and ABAQUS_DOWNLOAD_URL not set, unable to + download file. To enable automatic model downloading, set environment variable + ABAQUS_URL to point url to models.""") + return 1 + end + url = ENV["ABAQUS_DOWNLOAD_URL"] + if haskey(ENV, "ABAQUS_DOWNLOAD_DIR") + fn = rstrip(ENV["ABAQUS_DOWNLOAD_DIR"], '/') * "/" * fn + end + if !isfile(fn) + info("Downloading model $name from $url to $fn") + download("$url/$name.inp", fn) + end + return 0 +end + +""" Return input file name. """ +function abaqus_input_file_name(name) + fn = "$name.inp" + isfile(fn) && return fn + if haskey(ENV, "ABAQUS_DOWNLOAD_DIR") + fn = rstrip(ENV["ABAQUS_DOWNLOAD_DIR"], '/') * "/" * fn + end + isfile(fn) && return fn + return "" +end + +function abaqus_input_file_path(name) + return dirname(abaqus_input_file_name(name)) +end + +function abaqus_open_results(name) + path = abaqus_input_file_path(name) + result_file = "$path/$name.xmf" + return XDMF(result_file) end ### JuliaFEM-ABAQUS interface entry point -""" Run ABAQUS .inp file. """ -function abaqus_run_model(fn) - model = abaqus_read_model(fn) - model() -end +""" +Run ABAQUS model. If input file is not found, attempt to fetch it from internet +if fetch is set to true and ABAQUS_DOWNLOAD_URL is set. Return exit code 0 if +execution of model is success. +""" +function abaqus_run_model(name; fetch=false, verbose=false) -""" -Run ABAQUS .inp file. This function can be used to run some abaqus inp file -with the following extra feature: if file is not found from temporary directory, -attempt to download if from internet, if environment variable ABAQUS_TEST_URL is -set. This can be used to verify JuliaFEM code using well known ABAQUS test cases, -if they are found from company intranet and are accessible using wget / curl. -""" -function abaqus_run_test(name; print_test_file=false) - # get test file - fn = tempdir()*"/$name.inp" - if !isfile(fn) - # if file does not exist, attempt to download it if ABAQUS_TEST_URL is set - if haskey(ENV, "ABAQUS_TEST_URL") - url = ENV["ABAQUS_TEST_URL"]*"/$name.inp" - info("Downloading test from url $url") - download(url, fn) - else - # otherwise give up, no testing possible this time - info(""" - File $fn not found and ABAQUS_TEST_URL not set, unable to download file. - To access Abaqus test files, set environment variable ABAQUS_TEST_URL to - point url to Abaqus verification book. Typically the address is something - like `http://myurl.com:/books/eif`, so you need to set - `export ABAQUS_TEST_URL=http://myurl.com:/books/eif` to your .bashrc - file to make tests working.""") - return false - end + if !isfile("$name.inp") && fetch + status = abaqus_download(name) + status == 0 || return status # download failed end - if print_test_file + + fn = abaqus_input_file_name(name) + + if verbose + println(repeat("-", 80)) + println("Running ABAQUS model $name from file $fn") println(repeat("-", 80)) println(readall(fn)) println(repeat("-", 80)) end - abaqus_run_model(fn) + + model = abaqus_read_model(fn) + status = model() + return status end diff --git a/src/postprocess_utils.jl b/src/postprocess_utils.jl index f6e4539..61bf4f4 100644 --- a/src/postprocess_utils.jl +++ b/src/postprocess_utils.jl @@ -47,6 +47,31 @@ function convert(::Type{DataFrame}, dfs::AbstractString) return readtable(fn) end +function getindex(df::DataFrame, ids::Vector{Symbol}, cols::Vector{Symbol}) + rows = Int64[find(df[:id] .== id)[1] for id in ids] + return df[rows, cols] +end + +function getindex(df::DataFrame, id::Symbol, col::Symbol) + return getindex(df, [id], [col]) +end + +function getindex(df::DataFrame, id::Symbol, cols::Vector{Symbol}) + return getindex(df, [id], cols) +end + +function getindex(df::DataFrame, ids::Vector{Symbol}, col::Symbol) + return getindex(df, ids, [col]) +end + +function vec(df::DataFrame) + return vec(convert(Matrix{Float64}, df)) +end + +function isapprox(d1::DataFrame, d2::Vector) + return isapprox(vec(d1), d2) +end + """ Calculate field values to nodal points from Gauss points using least-squares fitting. """ @@ -186,6 +211,23 @@ function call(problem::Problem, field_name::AbstractString, time::Float64=0.0) return f end +function call(problem::Problem, ::Type{DataFrame}, field_name::AbstractString, + abbreviation::Symbol, time::Float64=0.0) + u = problem(field_name, time) + node_ids = collect(keys(u)) + column_names = [:id] + n = length(u[first(node_ids)]) + index = [Symbol("N$id") for id in node_ids] + result = Any[index] + for dof=1:n + push!(result, [u[id][dof] for id in node_ids]) + push!(column_names, Symbol("$abbreviation$dof")) + end + df = DataFrame(result, column_names) + sort!(df, cols=[:id]) + return df +end + """ Interpolate field from a set of elements. """ function call(problem::Problem, field_name::AbstractString, X::Vector, time::Float64=0.0; fillna=NaN) for element in get_elements(problem) diff --git a/src/postprocess_xdmf.jl b/src/postprocess_xdmf.jl index d76443e..c397524 100644 --- a/src/postprocess_xdmf.jl +++ b/src/postprocess_xdmf.jl @@ -1,8 +1,8 @@ # This file is a part of JuliaFEM. # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md -using LightXML using JuliaFEM +using LightXML # element codes: http://www.paraview.org/pipermail/paraview/2013-July/028859.html # > from ./VTK/ThirdParty/xdmf2/vtkxdmf2/libsrc/XdmfTopology.h diff --git a/src/problems.jl b/src/problems.jl index 0126356..c01d0f1 100644 --- a/src/problems.jl +++ b/src/problems.jl @@ -188,7 +188,7 @@ function initialize!(problem::Problem, time=0.0) end """ Update problem solution vector for assembly. """ -function update_assembly!(problem, u, la) +function update_assembly!(problem, u, la; verbose=false) assembly = get_assembly(problem) @@ -209,15 +209,15 @@ function update_assembly!(problem, u, la) assembly.u_prev = copy(assembly.u) assembly.la_prev = copy(assembly.la) if get_formulation_type(problem) == :total - info("$(problem.name): total formulation, replacing solution vector with new values") + verbose && info("$(problem.name): total formulation, replacing solution vector with new values") assembly.u = u assembly.la = la elseif get_formulation_type(problem) == :incremental - info("$(problem.name): incremental formulation, adding increment to solution vector") + verbose && info("$(problem.name): incremental formulation, adding increment to solution vector") assembly.u += u assembly.la = la elseif get_formulation_type(problem) == :forwarddiff - info("$(problem.name): forwarddiff formulation, adding increment to solution vector and reaction force vector") + verbose && info("$(problem.name): forwarddiff formulation, adding increment to solution vector and reaction force vector") assembly.u += u assembly.la += la else diff --git a/test/test_abaqus_verification.jl b/test/test_abaqus_verification.jl index 9a1a3c7..db7ac5a 100644 --- a/test/test_abaqus_verification.jl +++ b/test/test_abaqus_verification.jl @@ -7,6 +7,11 @@ using JuliaFEM.Postprocess using JuliaFEM.Abaqus using JuliaFEM.Testing +# to turn on automatic file download, set +# ENV["ABAQUS_DOWNLOAD_URL"] = http://:2080/v2016/books/eif +# if don't want to download all stuff to current directory, +# set also e.g. ENV["ABAQUS_DOWNLOAD_DIR"] = "/tmp" + #= test_name = "ecs4sfs1" @testset "$test_name" begin @@ -16,10 +21,12 @@ end =# @testset "ec38sfs2" begin - abaqus_run_test("ec38sfs2"; print_test_file=true) || return - #= - results = abaqus_read_results("ec38sfs2") - side = get_results(results, "SECTION"; name="side") + return_code = abaqus_run_model("ec38sfs2"; fetch=true, verbose=true) + return_code == 0 || return + @test return_code == 0 +#= + xdmf = abaqus_open_results("ec38sfs2") + side, opts = read_result(xdmf, "SECTION/side") @test isapprox(side["SOFM"], 3464.0) @test isapprox(side["SOF1"], 2000.0) @test isapprox(side["SOF2"], 2000.0) @@ -32,5 +39,5 @@ end @test isapprox(side["SOCF1"], 2/3) @test isapprox(side["SOCF2"], 2/3) @test isapprox(side["SOCF3"], 1/6) - =# +=# end diff --git a/test/test_elasticity_2d_linear_with_surface_load.jl b/test/test_elasticity_2d_linear_with_surface_load.jl index e19e7db..d582de3 100644 --- a/test/test_elasticity_2d_linear_with_surface_load.jl +++ b/test/test_elasticity_2d_linear_with_surface_load.jl @@ -5,7 +5,6 @@ using JuliaFEM using JuliaFEM.Preprocess using JuliaFEM.Postprocess using JuliaFEM.Testing -using JLD function JuliaFEM.get_model(::Type{Val{Symbol("test 2d linear elasticity with surface + volume load")}}) meshfile = "/geometry/2d_block/BLOCK_1elem.med" @@ -51,10 +50,15 @@ end E = 288.0 nu = 1/3 u3_expected = f/E*[-nu, 1] + g/(2*E)*[-nu, 1] - u3 = reshape(block.assembly.u, 2, 4)[:,3] - info("u3 = $u3") + + results = block(DataFrame, "displacement", :U, 0.0) + println(results) + u3 = results[:N3, [:U1, :U2]] + info("(u1,u2) at node 3") + info(u3) @test isapprox(u3, u3_expected) +#= info("strain") for ip in get_integration_points(block.elements[1]) eps = ip("strain") @@ -79,6 +83,8 @@ end # TODO: to postprocess #@test isapprox(stress[1], [0.0, g, 0.0]) #@test isapprox(strain[1], [u3[1], u3[2], 0.0]) +=# + end #= TODO: to other file diff --git a/test/test_postprocess.jl b/test/test_postprocess.jl index 703f37b..353210b 100644 --- a/test/test_postprocess.jl +++ b/test/test_postprocess.jl @@ -1,3 +1,6 @@ +# 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 diff --git a/test/test_xdmf.jl b/test/test_postprocess_xdmf.jl similarity index 100% rename from test/test_xdmf.jl rename to test/test_postprocess_xdmf.jl