improved postprocessing capabilities; use dataframes to show results

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. """

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:<port>/books/eif`, so you need to set
-            `export ABAQUS_TEST_URL=http://myurl.com:<port>/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
 

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)

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

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

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://<domain>: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

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

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