diff --git a/src/JuliaFEM.jl b/src/JuliaFEM.jl index b5c1353..7cca197 100644 --- a/src/JuliaFEM.jl +++ b/src/JuliaFEM.jl @@ -65,7 +65,8 @@ export AbstractSolver, Solver, Nonlinear, NonlinearSolver, Linear, LinearSolver, get_unknown_field_name, get_formulation_type, get_problems, get_field_problems, get_boundary_problems, get_field_assembly, get_boundary_assembly, - initialize!, create_projection, eliminate_interior_dofs + initialize!, create_projection, eliminate_interior_dofs, + is_field_problem, is_boundary_problem include("solvers_modal.jl") export Modal @@ -135,7 +136,8 @@ export calc_nodal_values!, copy_field!, calculate_area, calculate_center_of_mass, - calculate_second_moment_of_mass + calculate_second_moment_of_mass, + extract include("postprocess_xdmf.jl") export XDMF, xdmf_new_result!, xdmf_save_field!, xdmf_save! export DataFrame diff --git a/src/abaqus.jl b/src/abaqus.jl index b8294d3..a52dfbd 100644 --- a/src/abaqus.jl +++ b/src/abaqus.jl @@ -6,7 +6,6 @@ importall Base using JuliaFEM using JuliaFEM.Preprocess using JuliaFEM.Postprocess -using LightXML ### Model definitions for ABAQUS data model @@ -265,8 +264,9 @@ end typealias BOUNDARY_CONDITIONS Union{BOUNDARY, CLOAD, DLOAD, DSLOAD} @register_abaqus_keyword("NODE PRINT") +@register_abaqus_keyword("EL PRINT") @register_abaqus_keyword("SECTION PRINT") -typealias OUTPUT_REQUESTS Union{NODE_PRINT, SECTION_PRINT} +typealias OUTPUT_REQUESTS Union{NODE_PRINT, EL_PRINT, SECTION_PRINT} ## Properties @@ -391,7 +391,7 @@ end """ Dirichlet boundary condition. """ function create_boundary_problem(model::Model, bc::AbstractBoundaryCondition, ::BOUNDARY; verbose=true) dim = determine_problem_dimension(model) - problem = Problem(Dirichlet, "Dirichlet bc *BOUNDARY", dim, "displacement") + problem = Problem(Dirichlet, "Dirichlet boundary *BOUNDARY", dim, "displacement") for row in bc.data if isa(row[1], AbstractString) # node set given @@ -462,7 +462,7 @@ function create_boundary_problem(model::Model, bc::AbstractBoundaryCondition, :: child_element = Element(JuliaFEM.(child_element_type), child_element_connectivity) update!(child_element, "geometry", model.mesh.nodes) - update!(child_element, "surface pressure", pressure) + update!(child_element, "surface pressure", -pressure) push!(problem.elements, child_element) end return problem @@ -472,13 +472,18 @@ end function create_boundary_problem(model::Model, bc::AbstractBoundaryCondition, ::CLOAD; verbose=false) dim = determine_problem_dimension(model) problem = Problem(Elasticity, "Concentrated load *CLOAD", dim) - for row in bc.data - node, dof, load = row + nodes = sort(unique([row[1] for row in bc.data])) + elements = Dict() + for node in nodes element = Element(Poi1, [node]) update!(element, "geometry", model.mesh.nodes) - update!(element, "displacement traction force $dof", load) - push!(problem.elements, element) + elements[node] = element end + for row in bc.data + node, dof, load = row + update!(elements[node], "concentrated force $dof", load) + end + problem.elements = collect(values(elements)) return problem end @@ -531,8 +536,6 @@ function get_child_element(element_type::Symbol, element_side::Symbol, process_output_request(model, solver, output_request, kind, target) end -using DataFrames - function process_output_request(model::Model, solver::Solver, output_request::AbstractOutputRequest, ::Type{Val{:NODE}}, ::Type{Val{:PRINT}}) data = output_request.data @@ -546,7 +549,7 @@ function process_output_request(model::Model, solver::Solver, output_request::Ab for row in data info(repeat("-", 80)) codes = join(row, ", ") - info("*NODE OUTPUT request, with fields $codes") + info("*NODE PRINT request, with fields $codes") if length(options) != 0 info("Additional options: $options") end @@ -554,15 +557,56 @@ function process_output_request(model::Model, solver::Solver, output_request::Ab tables = Any[] for code in row haskey(code_mapping, code) || continue - for problem in model.problems - field_name = code_mapping[code] - abbr = get(abbr_mapping, code, code) - table = problem(DataFrame, field_name, abbr, solver.time) - push!(tables, table) - end + field_name = code_mapping[code] + abbr = get(abbr_mapping, code, code) + table = solver(DataFrame, field_name, abbr, solver.time) + push!(tables, table) end length(tables) != 0 || continue - results = join(tables..., on=:id, kind=:outer) + results = join(tables..., on=:NODE, kind=:outer) + sort!(results, cols=[:NODE]) + println() + println(results) + println() + end +end + +function process_output_request(model::Model, solver::Solver, output_request::AbstractOutputRequest, + ::Type{Val{:EL}}, ::Type{Val{:PRINT}}) + data = output_request.data + options = output_request.options + code_mapping = Dict( + :COORD => "geometry", + :S => "stress", + :E => "strain") + abbr_mapping = Dict(:COORD => :COOR) + for row in data + info(repeat("-", 80)) + codes = join(row, ", ") + info("*EL PRINT request, with fields $codes") + if length(options) != 0 + info("Additional options: $options") + end + info(repeat("-", 80)) + tables = Any[] + for code in row + haskey(code_mapping, code) || continue + field_name = code_mapping[code] + abbr = get(abbr_mapping, code, code) + table = solver(DataFrame, solver.time, Val{code}) + push!(tables, table) + end + length(tables) != 0 || continue + results = first(tables) + if length(tables) > 1 + for i=2:length(tables) + results = join(results, tables[i], on=:ELEMENT, kind=:outer) + end + end + sort!(results, cols=[:ELEMENT, :IP]) + # filter out elements with id -1, they are automatically created boundary elements + fel = find(results[:ELEMENT] .!= Symbol("E-1")) + results = results[fel, :] println() println(results) println() diff --git a/src/elements.jl b/src/elements.jl index ee72ea4..17c46dd 100644 --- a/src/elements.jl +++ b/src/elements.jl @@ -168,11 +168,13 @@ function update!(element::Element, field_name::AbstractString, datas::Union{Real end end -function update!(element::Element, field_name, datas::Pair...) +#= +function update!(element::Element, field_name, data::Pair...) for data in datas update!(element, field_name, data) end end +=# function update!(element::Element, field_name, data::Pair{Float64, Vector{Any}}) if haskey(element, field_name) @@ -190,6 +192,14 @@ function update!(element::Element, field_name, data::Pair{Float64, Vector{Int64} end end +function update!(element::Element, field_name, data::Pair{Float64, Vector{Float64}}) + if haskey(element, field_name) + update!(element[field_name], data) + else + element[field_name] = data + end +end + function update!(element::Element, field_name, data::Pair{Float64, Vector{Vector{Float64}}}) if haskey(element, field_name) update!(element[field_name], data) diff --git a/src/elements_lagrange.jl b/src/elements_lagrange.jl index a788c3b..622a1b4 100644 --- a/src/elements_lagrange.jl +++ b/src/elements_lagrange.jl @@ -22,6 +22,10 @@ function get_basis(element::Element{Poi1}, ip, time) return [1] end +function get_dbasis(element::Element{Poi1}, ip, time) + return [0] +end + function call(element::Element{Poi1}, ip, time::Float64, ::Type{Val{:detJ}}) return 1.0 end diff --git a/src/postprocess_utils.jl b/src/postprocess_utils.jl index 61bf4f4..d5c0eef 100644 --- a/src/postprocess_utils.jl +++ b/src/postprocess_utils.jl @@ -7,6 +7,7 @@ using JuliaFEM using DataFrames using HDF5 using LightXML +using StringUtils import HDF5: h5read, h5write @@ -47,21 +48,14 @@ 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]) +function extract(df::DataFrame, args...; kwargs...) + result = copy(df) + for (k,v) in kwargs + rows = find(df[k] .== v) + result = result[rows, :] + end + foo = Symbol[si for si in args] + return result[foo] end function vec(df::DataFrame) @@ -72,6 +66,12 @@ function isapprox(d1::DataFrame, d2::Vector) return isapprox(vec(d1), d2) end +""" A more appropriate representation for floats in results. """ +function DataFrames.ourshowcompact(io::IO, x::Float64) + print(io, u"\% 0.4E(x)") + return +end + """ Calculate field values to nodal points from Gauss points using least-squares fitting. """ @@ -201,9 +201,13 @@ end function call(problem::Problem, field_name::AbstractString, time::Float64=0.0) f = Dict() for element in get_elements(problem) + haskey(element, field_name) || continue for (c, v) in zip(get_connectivity(element), element(field_name, time)) if haskey(f, c) - @assert isapprox(f[c], v) + 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 end f[c] = v end @@ -211,11 +215,10 @@ 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) +function to_dataframe(u::Dict, abbreviation::Symbol) + length(u) != 0 || return DataFrame() node_ids = collect(keys(u)) - column_names = [:id] + column_names = [:NODE] n = length(u[first(node_ids)]) index = [Symbol("N$id") for id in node_ids] result = Any[index] @@ -224,10 +227,92 @@ function call(problem::Problem, ::Type{DataFrame}, field_name::AbstractString, push!(column_names, Symbol("$abbreviation$dof")) end df = DataFrame(result, column_names) - sort!(df, cols=[:id]) + sort!(df, cols=[:NODE]) return df end +function call(problem::Problem, ::Type{DataFrame}, field_name::AbstractString, + abbreviation::Symbol, time::Float64=0.0) + u = problem(field_name, time) + return to_dataframe(u, abbreviation) +end + +function call(solver::Solver, ::Type{DataFrame}, field_name::AbstractString, + abbreviation::Symbol, time::Float64=0.0) + u = Dict() + for problem in get_problems(solver) + u = merge(u, problem(field_name, time)) + end + return to_dataframe(u, abbreviation) +end + +function get_components(n, m) + if n == m + if n == 1 + return Vector{Int}[[1,1]] + end + if n == 2 + return Vector{Int}[[1,1], [2,2], [1,2]] + elseif n == 3 + return Vector{Int}[[1,1], [2,2], [3,3], [1,2], [1,3], [2,3]] + else + error("get_components, n=$n, m=$m!") + end + end +end + +""" Return T in integration points. """ +function call{T}(problem::Problem, ::Type{DataFrame}, element::Element, time::Float64, + ::Type{Val{T}}) + column_names = [:ELEMENT, :IP] + ips = get_integration_points(element) + field = Any[problem(element, ip, time, Val{T}) for ip in ips] + # FIXME, handle better ..? + first(field) == nothing && return DataFrame() + m = length(field) + n = length(first(field)) + result = Any[] + push!(result, [Symbol("E$(element.id)") for i=1:m]) + push!(result, [Symbol("P$i") for i=1:m]) + is_tensor_field = isa(first(field), Matrix) + if is_tensor_field + n, m = size(first(field)) + components = get_components(n, m) + for (j, k) in components + push!(column_names, Symbol("$T$j$k")) + push!(result, [S[j,k] for S in field]) + end + else + components = collect(1:n) + for j in components + push!(column_names, Symbol("$T$j")) + push!(result, [S[j] for S in field]) + end + end + df = DataFrame(result, column_names) + sort!(df, cols=[:IP]) +end + +function call{T}(problem::Problem, ::Type{DataFrame}, time::Float64, ::Type{Val{T}}) + tables = [problem(DataFrame, element, time, Val{T}) for element in get_elements(problem)] + results = [tables...;] + return results +end + +function call{T}(solver::Solver, ::Type{DataFrame}, time::Float64, ::Type{Val{T}}) + problems = get_problems(solver) + tables = Any[] + for problem in get_problems(solver) + try + push!(tables, problem(DataFrame, time, Val{T})) + catch + warn("Unable to obtain results $T for problem $(problem.name)") + end + end + results = [tables...;] + return results +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/preprocess.jl b/src/preprocess.jl index bda6cc7..95a0375 100644 --- a/src/preprocess.jl +++ b/src/preprocess.jl @@ -91,13 +91,21 @@ function filter_by_element_set(mesh::Mesh, set_name) filter_by_element_id(mesh::Mesh, collect(mesh.element_sets[set_name])) end +function create_element(mesh::Mesh, id::Int) + connectivity = mesh.elements[id] + element_type = JuliaFEM.(mesh.element_types[id]) + element = Element(element_type, connectivity) + update!(element, "geometry", mesh.nodes) + element.id = id + return element +end + function create_elements(mesh::Mesh; element_type=nothing) element_ids = collect(keys(mesh.elements)) if element_type != nothing filter!(id -> mesh.element_types[id] == element_type, element_ids) end - elements = [Element(JuliaFEM.(mesh.element_types[id]), mesh.elements[id]) for id in element_ids] - update!(elements, "geometry", mesh.nodes) + elements = [create_element(mesh, id) for id in element_ids] return elements end @@ -115,8 +123,7 @@ function create_elements(mesh::Mesh, element_sets::Symbol...; element_type=nothi filter!(id -> mesh.element_types[id] == element_type, element_ids) end - elements = [Element(JuliaFEM.(mesh.element_types[id]), mesh.elements[id]) for id in element_ids] - update!(elements, "geometry", mesh.nodes) + elements = [create_element(mesh, id) for id in element_ids] return elements end diff --git a/src/problems.jl b/src/problems.jl index c01d0f1..0e36e7a 100644 --- a/src/problems.jl +++ b/src/problems.jl @@ -141,63 +141,53 @@ function get_assembly(problem) return problem.assembly end -""" Initialize unknown field ready for nonlinear iterations, i.e., - take last known value and set it as a initial quess for next - time increment. -""" -function initialize!(problem::Problem, time=0.0) + +""" Initialize element ready for calculation. """ +function initialize!(problem::Problem, element::Element, time::Float64) field_name = get_unknown_field_name(problem) field_dim = get_unknown_field_dimension(problem) - for element in get_elements(problem) - gdofs = get_gdofs(problem, element) - if haskey(element, field_name) - # if field is found, copy last known solution to new time as initial guess - field = last(element[field_name]) - if !isa(field, TimeVariantField) - info("Unable to initialize field $field_name for problem, is not time variant?") - continue - end + nnodes = length(element) - if !isapprox(field.time, time) - last_data = copy(last(element[field_name]).data) - push!(element[field_name], time => last_data) - end - else # if field not found at all, initialize new zero field. - data = Vector{Float64}[zeros(field_dim) for i in 1:length(element)] - element[field_name] = (time => data) + # initialize primary field + if !haskey(element, field_name) + if field_dim == 1 + update!(element, field_name, time => zeros(nnodes)) + else + update!(element, field_name, time => [zeros(field_dim) for i=1:nnodes]) end end - # if this is boundary problem and not dirichlet problem, initialize field - # for primary variable too + + # if boundary problem, initialize field for main problem too is_boundary_problem(problem) || return - #is_dirichlet_problem(problem) && return field_name = get_parent_field_name(problem) - for element in get_elements(problem) - gdofs = get_gdofs(problem, element) - if haskey(element, field_name) - # if field is found, copy last known solution to new time as initial guess - if !isapprox(last(element[field_name]).time, time) - last_data = copy(last(element[field_name]).data) - push!(element[field_name], time => last_data) - end - else # if field not found at all, initialize new zero field. - data = Vector{Float64}[zeros(field_dim) for i in 1:length(element)] - element[field_name] = (time => data) + if !haskey(element, field_name) + if field_dim == 1 + update!(element, field_name, time => zeros(nnodes)) + else + update!(element, field_name, time => [zeros(field_dim) for i=1:nnodes]) end end end -""" Update problem solution vector for assembly. """ -function update_assembly!(problem, u, la; verbose=false) +function initialize!(problem::Problem, time::Float64=0.0) + for element in get_elements(problem) + initialize!(problem, element, time) + end +end - assembly = get_assembly(problem) +""" Update problem solution vector for assembly. """ +function update!(problem::Problem, assembly::Assembly, u::Vector, la::Vector; verbose=false) # resize & fill with zeros vectors if length mismatch with current solution + if length(u) != length(assembly.u) + info("resizing solution vector u") resize!(assembly.u, length(u)) fill!(assembly.u, 0.0) end + if length(la) != length(assembly.la) + info("resizing lagrange multipliers vector u") resize!(assembly.la, length(la)) fill!(assembly.la, 0.0) end @@ -228,47 +218,54 @@ function update_assembly!(problem, u, la; verbose=false) # calculate change of norm assembly.u_norm_change = norm(assembly.u - assembly.u_prev) assembly.la_norm_change = norm(assembly.la - assembly.la_prev) - #return assembly.u_norm_change, assembly.la_norm_change return assembly.u, assembly.la end -""" Update solutions to elements. +""" Return global solution (u, la) for problem. Notes ----- -This assumes that element is properly initialized so that last known field data -is from current time. For boundary problems solution is updated from lambda vector -and for field problems from actual solution vector. +If length of solution vector != number of nodes, i.e. field dimension is +something other than 1, reshape vectors so it's length matches to the +number of nodes so that one can easily get nodal results. """ -function update_elements!{P<:FieldProblem}(problem::Problem{P}, u, la) - field_name = get_unknown_field_name(problem) +function get_global_solution(problem::Problem, assembly::Assembly) + u = assembly.u + la = assembly.la field_dim = get_unknown_field_dimension(problem) - nnodes = round(Int, length(u)/field_dim) - solution = reshape(u, field_dim, nnodes) - for element in get_elements(problem) - connectivity = get_connectivity(element) # node ids - local_sol = Vector{Float64}[solution[:, node_id] for node_id in connectivity] - last(element[field_name]).data = local_sol + if field_dim == 1 + return u, la + else + nnodes = round(Int, length(u)/field_dim) + u = reshape(u, field_dim, nnodes) + u = Vector{Float64}[u[:,i] for i in 1:nnodes] + la = reshape(la, field_dim, nnodes) + la = Vector{Float64}[la[:,i] for i in 1:nnodes] + return u, la end end -function update_elements!{P<:BoundaryProblem}(problem::Problem{P}, u, la) + +""" Update solution from assebly to elements. """ +function update!{P<:FieldProblem}(problem::Problem{P}, assembly::Assembly, elements::Vector{Element}, time::Float64) + u, la = get_global_solution(problem, assembly) field_name = get_unknown_field_name(problem) - field_dim = get_unknown_field_dimension(problem) - nnodes = round(Int, length(u)/field_dim) - solution = reshape(la, field_dim, nnodes) - for element in get_elements(problem) - connectivity = get_connectivity(element) # node ids - local_sol = Vector{Float64}[solution[:, node_id] for node_id in connectivity] - last(element[field_name]).data = local_sol + # update solution u for elements + for element in elements + connectivity = get_connectivity(element) + update!(element, field_name, time => u[connectivity]) end - # if boundary problem is not dirichlet, update also data of main problem - # is_dirichlet_problem(problem) && return - field_name = get_parent_field_name(problem) - solution = reshape(u, field_dim, nnodes) - for element in get_elements(problem) - connectivity = get_connectivity(element) # node ids - local_sol = Vector{Float64}[solution[:, node_id] for node_id in connectivity] - last(element[field_name]).data = local_sol +end + +function update!{P<:BoundaryProblem}(problem::Problem{P}, assembly::Assembly, elements::Vector{Element}, time::Float64) + u, la = get_global_solution(problem, assembly) + parent_field_name = get_parent_field_name(problem) # displacement + field_name = get_unknown_field_name(problem) # reaction force + # update solution u and reaction force λ for boundary elements + for element in elements + connectivity = get_connectivity(element) + update!(element, parent_field_name, time => u[connectivity]) + # FIXME + update!(element, field_name, time => -la[connectivity]) end end @@ -276,12 +273,16 @@ function get_elements(problem::Problem) return problem.elements end +function get_assembly(problem::Problem) + return problem.assembly +end + function length(problem::Problem) return length(problem.elements) end -function update!(problem::Problem, field_name, field) - update!(problem.elements, field_name, field) +function update!(problem::Problem, field_name::AbstractString, data) + update!(problem.elements, field_name::AbstractString, data) end """ Return the dimension of the unknown field of this problem. """ @@ -380,3 +381,4 @@ function find_nodes_by_dofs(dim, dofs) end return nodes end + diff --git a/src/problems_elasticity.jl b/src/problems_elasticity.jl index d8ff22a..81c5a71 100644 --- a/src/problems_elasticity.jl +++ b/src/problems_elasticity.jl @@ -530,6 +530,7 @@ function assemble{El<:Elasticity3DSurfaceElements}(problem::Problem{Elasticity}, Kg = zeros(dim*nnodes, dim*nnodes) f = zeros(dim*nnodes) + has_concentrated_forces = false for ip in get_integration_points(element) detJ = element(ip, time, Val{:detJ}) w = ip.weight*detJ @@ -543,6 +544,11 @@ function assemble{El<:Elasticity3DSurfaceElements}(problem::Problem{Elasticity}, T = element("displacement traction force $i", ip, time) f[i:dim:end] += w*vec(T*N) end + if haskey(element, "concentrated force $i") + has_concentrated_forces = true + T = element("concentrated force $i", ip, time) + f[i:dim:end] += w*vec(T*N) + end end if haskey(element, "surface pressure") J = element(ip, time, Val{:Jacobian})' @@ -553,6 +559,9 @@ function assemble{El<:Elasticity3DSurfaceElements}(problem::Problem{Elasticity}, f += w*p*vec(n*N) end end + if has_concentrated_forces + update!(element, "concentrated force", time => Any[f]) + end return Km, Kg, f end @@ -710,3 +719,32 @@ end =# + +function call(problem::Problem, element::Element, ip, time::Float64, ::Type{Val{:E}}) + haskey(element, "displacement") || return nothing + gradu = element("displacement", ip, time, Val{:Grad}) + eps = 0.5*(gradu + gradu') + return eps +end + +function call(problem::Problem, element::Element, ip, time::Float64, ::Type{Val{:S}}) + haskey(element, "displacement") || return nothing + props = problem.properties + eps = problem(element, ip, time, Val{:E}) + eps == nothing && return nothing + E = element("youngs modulus", ip, time) + nu = element("poissons ratio", ip, time) + mu = E/(2.0*(1.0+nu)) + la = E*nu/((1.0+nu)*(1.0-2.0*nu)) + if props.formulation in [:plane_stress, :plane_strain] + la = 2.0*la*mu/(la+2.0*mu) + end + S = la*trace(eps)*I + 2.0*mu*eps + return S +end + +function call(problem::Problem, element::Element, ip, time::Float64, ::Type{Val{:COORD}}) + haskey(element, "geometry") || return nothing + return element("geometry", ip, time) +end + diff --git a/src/solvers.jl b/src/solvers.jl index 93b8c48..fe7da6d 100644 --- a/src/solvers.jl +++ b/src/solvers.jl @@ -408,8 +408,12 @@ function update!(solver::Solver, u::Vector, la::Vector; show_info=true) show_info && info("Updating problems ...") t0 = Base.time() for problem in solver.problems - u_new, la_new = update_assembly!(problem, u, la) - update_elements!(problem, u_new, la_new) + assembly = get_assembly(problem) + elements = get_elements(problem) + # update solution, first for assembly (u,la) ... + update!(problem, assembly, u, la) + # .. and then from assembly (u,la) to elements + update!(problem, assembly, elements, solver.time) end t1 = round(Base.time()-t0, 2) show_info && info("Updated problems in $t1 seconds.") diff --git a/src/solvers_modal.jl b/src/solvers_modal.jl index 933a4b6..a298c2e 100644 --- a/src/solvers_modal.jl +++ b/src/solvers_modal.jl @@ -101,21 +101,23 @@ function call(solver::Solver{Modal}; show_info=true, debug=false) info("Eigenvalues computed in $t1 seconds. Eigenvalues: $om2") for i=1:length(om2) + freq = real(sqrt(om2[i])/(2.0*pi)) u = props.eigvecs[:,i] field_dim = get_unknown_field_dimension(solver) field_name = get_unknown_field_name(solver) - nnodes = round(Int, length(u)/field_dim) - solution = reshape(u, field_dim, nnodes) + if field_dim != 1 + nnodes = round(Int, length(u)/field_dim) + u = reshape(u, field_dim, nnodes) + u = Vector{Float64}[u[:,i] for i in 1:nnodes] + end for problem in get_problems(solver) - local_sol = Dict{Int64, Vector{Float64}}() - for node_id in get_connectivity(problem) - local_sol[node_id] = solution[:, node_id] + for element in get_elements(problem) + connectivity = get_connectivity(element) + update!(element, field_name, freq => u[connectivity]) end - freq = real(sqrt(om2[i])/(2.0*pi)) - update!(problem, field_name, freq => local_sol) end end - + return true end diff --git a/test/test_abaqus_verification.jl b/test/test_abaqus_verification.jl index db7ac5a..5225af1 100644 --- a/test/test_abaqus_verification.jl +++ b/test/test_abaqus_verification.jl @@ -8,36 +8,64 @@ 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" +# ENV["ABAQUS_DOWNLOAD_URL"] = "http://:2080/v2016/books/eif" +# if don't want to download all stuff to current directory, set also +# ENV["ABAQUS_DOWNLOAD_DIR"] = "/tmp" -#= -test_name = "ecs4sfs1" -@testset "$test_name" begin - abaqus_run_test(test_name) || return - results = abaqus_read_results(test_name) +""" Run test, return true if simulation is succesfull, i.e. no errors raise +during parsing .inp file or execution of model. This doesn't mean that results +are meaningful; they must be checked in separately. Running model only verifies +that no catastrophic failures happen during file parsing. """ +function abaqus_run_test(name) + return_code = abaqus_run_model(name; fetch=true, verbose=true) + return_code == 0 && return true + return false end -=# -@testset "ec38sfs2" begin - 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) - @test isapprox(side["SOF3"], 2000.0) - @test isapprox(side["SOMM"], 2828.0) - @test isapprox(side["SOM1"], 0.0) - @test isapprox(side["SOM2"], 2000.0) - @test isapprox(side["SOM3"], -2000.0) - @test isapprox(side["SOAREA"], 2.000) - @test isapprox(side["SOCF1"], 2/3) - @test isapprox(side["SOCF2"], 2/3) - @test isapprox(side["SOCF3"], 1/6) -=# +@testset "JuliaFEM-ABAQUS interface" begin + @testset "1 Element Verification" begin + @testset "1.2 Eigenvalue tests" begin + @testset "1.2.1 Eigenvalue extraction for single unconstrained elements" begin + @testset "Acoustic elements" begin + @testset "AC1D2 elements." begin +# abaqus_run_test("ec12afe1") || return + end + end + @testset "Three-dimensional continuum elements" begin + @testset "C3D10 elements." begin +# abaqus_run_test("ec3asfe1") || return + end + end + end + end + @testset "1.3 Simple load tests" begin + @testset "1.3.1 Membrane loading of plane stress, plane strain, membrane, and shell elements" begin + @testset "CPS4 elements." begin +# abaqus_run_test("ecs4sfs1") || return + end + end + @testset "1.3.3 Three-dimensional solid elements" begin + @testset "C3D8 elements." begin + abaqus_run_test("ec38sfs2") || return + #= to check also results: + 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) + @test isapprox(side["SOF3"], 2000.0) + @test isapprox(side["SOMM"], 2828.0) + @test isapprox(side["SOM1"], 0.0) + @test isapprox(side["SOM2"], 2000.0) + @test isapprox(side["SOM3"], -2000.0) + @test isapprox(side["SOAREA"], 2.000) + @test isapprox(side["SOCF1"], 2/3) + @test isapprox(side["SOCF2"], 2/3) + @test isapprox(side["SOCF3"], 1/6) + =# + end + end + end + end end + diff --git a/test/test_contact_2d_small_sliding.jl b/test/test_contact_2d_small_sliding.jl index 7848520..484bd8e 100644 --- a/test/test_contact_2d_small_sliding.jl +++ b/test/test_contact_2d_small_sliding.jl @@ -122,7 +122,7 @@ end slaves = get_slave_elements(contact) node_ids, la = get_nodal_vector(slaves, "reaction force", 0.0) node_ids, n = get_nodal_vector(slaves, "normal", 0.0) - pres = [dot(ni, lai) for (ni, lai) in zip(n, la)] + pres = [dot(ni, -lai) for (ni, lai) in zip(n, la)] #@test isapprox(maximum(pres), 4060.010799583303) # 12 % error in maximum pressure @test isapprox(maximum(pres), 3585.0; rtol = 12.0e-2) @@ -137,7 +137,7 @@ end n = sel("normal", ip, time) t = Q'*n la = sel("reaction force", ip, time) - Rn += w*dot(n, la) + Rn += w*dot(n, -la) Rt += w*dot(t, la) end end diff --git a/test/test_elasticity_1d.jl b/test/test_elasticity_1d.jl new file mode 100644 index 0000000..06f61d0 --- /dev/null +++ b/test/test_elasticity_1d.jl @@ -0,0 +1,26 @@ +# 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 + +@testset "1d strain" begin + X = Dict{Int64, Vector{Float64}}( + 1 => [0.0, 0.0, 0.0], + 2 => [1.0, 1.0, 1.0]) + u = Dict{Int64, Vector{Float64}}( + 1 => [0.0, 0.0, 0.0], + 2 => [1.0, 1.0, 1.0]) + element = Element(Seg2, [1, 2]) + update!(element, "geometry", X) + detJ = element([0.0], 0.0, Val{:detJ}) + info("detJ = $detJ") + @test isapprox(detJ, sqrt(3)/2) + J = element([0.0], 0.0, Val{:Jacobian}) + info("J = $J") + @test isapprox(J, [0.5 0.5 0.5]) + update!(element, "displacement", u) + # FIXME +# gradu = element("displacement", [0.0], 0.0, Val{:Grad}) +# info("1d bar: ∇u = $gradu") +end diff --git a/test/test_elasticity_2d_linear_with_surface_load.jl b/test/test_elasticity_2d_linear_with_surface_load.jl index 7a2e9b5..6271663 100644 --- a/test/test_elasticity_2d_linear_with_surface_load.jl +++ b/test/test_elasticity_2d_linear_with_surface_load.jl @@ -7,9 +7,15 @@ using JuliaFEM.Postprocess using JuliaFEM.Testing #= -- solve 2d plane stress problem with known solution -- test postprocessing of nodal fields: (geometry, displacement - reaction force, concentrated force) +- solve 2d plane stress problem with known solution: + surface traction force in 2d + volume load in 2d + reaction force +- test postprocessing of nodal fields: + geometry + displacement + reaction force + concentrated force =# @testset "test 2d linear elasticity with surface + volume load" begin meshfile = "/geometry/2d_block/BLOCK_1elem.med" @@ -42,8 +48,13 @@ using JuliaFEM.Testing update!(bc_sym_13, "displacement 2", 0.0) solver = LinearSolver(block, traction, bc_sym_23, bc_sym_13) +# assemble!(solver) +# dump(full(bc_sym_23.assembly.C1)) solver() + info("u = ", block.assembly.u) + info("λ = ", block.assembly.la) + f = 288.0 g = 576.0 E = 288.0 @@ -51,40 +62,41 @@ using JuliaFEM.Testing u3_expected = f/E*[-nu, 1] + g/(2*E)*[-nu, 1] # fetch nodal results X + u and join them into one table using DataFrames - X = block(DataFrame, "geometry", :COOR, 0.0) - u = block(DataFrame, "displacement", :U, 0.0) - results = join(X, u, on=:id, kind=:outer) + X = solver(DataFrame, "geometry", :COOR) + u = solver(DataFrame, "displacement", :U) + la = solver(DataFrame, "reaction force", :RF) + f = solver(DataFrame, "concentrated force", :CF) + results = join(X, u, on=:NODE, kind=:outer) + results = join(results, la, on=:NODE, kind=:outer) + length(f) != 0 && (results = join(results, f, on=:NODE, kind=:outer)) + sort!(results, cols=[:NODE]) println(results) - u3 = results[:N3, [:U1, :U2]] + u3 = extract(results, NODE=:N3, :U1, :U2) @test isapprox(u3, u3_expected) -#= - info("strain") - for ip in get_integration_points(block.elements[1]) - eps = ip("strain") - @printf "%i | %8.3f %8.3f | %8.3f %8.3f %8.3f\n" ip.id ip.coords[1] ip.coords[2] eps[1] eps[2] eps[3] - # TODO: to postprocess ...? - #@test isapprox(eps, [u3[1], u3[2], 0.0]) - end + # element details + el = first(block.elements) + S1 = block(el, [0.0, 0.0], 0.0, Val{:S}) + S1 = S1[[1,4,2]] + E1= block(el, [0.0, 0.0], 0.0, Val{:E}) + E1 = E1[[1,4,2]] + C1 = block(el, [0.0, 0.0], 0.0, Val{:COORD}) + info("strain = $E1, stress = $S1, at $C1") + @test isapprox(E1, [-2/3, 2.0, 0.0]) + @test isapprox(S1, [0.0, 576.0, 0.0]) + @test isapprox(C1, [0.5, 0.5]) - info("stress") - for ip in get_integration_points(block.elements[1]) - sig = ip("stress") - @printf "%i | %8.3f %8.3f | %8.3f %8.3f %8.3f\n" ip.id ip.coords[1] ip.coords[2] sig[1] sig[2] sig[3] - # TODO: to postprocess - #@test isapprox(sig, [0.0, g, 0.0]) - end + S1 = block(DataFrame, 0.0, Val{:S}) + E1 = block(DataFrame, 0.0, Val{:E}) + C1 = block(DataFrame, 0.0, Val{:COORD}) - calc_nodal_values!(block.elements, "strain", 3, 0.0) - calc_nodal_values!(block.elements, "stress", 3, 0.0) - info(block.elements[1]["stress"](0.0)) - node_ids, strain = get_nodal_vector(block.elements, "strain", 0.0) - node_ids, stress = get_nodal_vector(block.elements, "stress", 0.0) - # TODO: to postprocess - #@test isapprox(stress[1], [0.0, g, 0.0]) - #@test isapprox(strain[1], [u3[1], u3[2], 0.0]) -=# + println(S1) + println(E1) + println(C1) + + S = solver(DataFrame, 0.0, Val{:S}) + println(S) end diff --git a/test/test_elements.jl b/test/test_elements.jl index 3017934..f2e10fa 100644 --- a/test/test_elements.jl +++ b/test/test_elements.jl @@ -106,6 +106,7 @@ end @test isapprox(fb, 1.0) end +#= unnecessary feature @testset "add two time dependent fields to element at once" begin el = Element(Seg2, [1, 2]) update!(el, "foo1", 1.0 => 1.0) @@ -113,6 +114,7 @@ end update!(el, "foo2", 1.0 => 1.0, 2.0 => 2.0) @test isapprox(el("foo1", 1.5), el("foo2", 1.5)) end +=# @testset "add elements to elements" begin el1 = Element(Seg2, [1, 2]) diff --git a/test/test_heat.jl b/test/test_heat.jl index 652ac09..c9bb627 100644 --- a/test/test_heat.jl +++ b/test/test_heat.jl @@ -58,7 +58,8 @@ end el2 = Element(Seg2, [1, 2]) update!(el2, "geometry", X) # linear ramp from 0 -> 6 in time 0 -> 1 - update!(el2, "temperature flux", 0.0 => 0.0, 1.0 => 6.0) + update!(el2, "temperature flux", 0.0 => 0.0) + update!(el2, "temperature flux", 1.0 => 6.0) # define heat problem and push elements to problem problem = Problem(Heat, "one element heat problem", 1) diff --git a/test/test_modal_analysis.jl b/test/test_modal_analysis.jl index 5654b3b..8f297de 100644 --- a/test/test_modal_analysis.jl +++ b/test/test_modal_analysis.jl @@ -4,7 +4,7 @@ using JuliaFEM using JuliaFEM.Testing -@testset "test eigenvalues for single tet4 element" begin +function get_model() X = Dict{Int, Vector{Float64}}( 1 => [2.0, 3.0, 4.0], 2 => [6.0, 3.0, 2.0], @@ -19,33 +19,38 @@ using JuliaFEM.Testing e2 = Element(Tri3, [1, 2, 3]) update!([e1, e2], "geometry", X) update!([e1, e2], "displacement", 0.0 => u) - update!(e1, "youngs modulus" => 96.0, - "poissons ratio" => 1.0/3.0, - "density" => 420.0) - update!(e2, "displacement 1" => 0.0, - "displacement 2" => 0.0, - "displacement 3" => 0.0) + update!(e1, "youngs modulus" => 96.0) + update!(e1, "poissons ratio" => 1.0/3.0) + update!(e1, "density" => 420.0) + update!(e2, "displacement 1" => 0.0) + update!(e2, "displacement 2" => 0.0) + update!(e2, "displacement 3" => 0.0) p1 = Problem(Elasticity, 3) p1.properties.finite_strain = false p1.properties.geometric_stiffness = false p2 = Problem(Dirichlet, p1) push!(p1, e1) push!(p2, e2) - s1 = Solver(Modal) - s1.properties.which = :LM - push!(s1, p1, p2) + solver = Solver(Modal) + solver.properties.which = :LM + push!(solver, p1, p2) + return solver +end - s1(; debug=true) - @test isapprox(s1.properties.eigvals, [4/3, 1/3]) +@testset "test eigenvalues for single tet4 element" begin + solver = get_model() + solver(; debug=true) + @test isapprox(solver.properties.eigvals, [4/3, 1/3]) +end - empty!(p1) - empty!(p2) - empty!(p1.assembly.M) -# p1.properties.finite_strain = true - p1.properties.geometric_stiffness = true - s1.properties.geometric_stiffness = true - s1(; debug=true) - @test isapprox(s1.properties.eigvals, [5/3, 2/3]) +@testset "test eigenvalues for single tet4 element, with geometric stiffness" begin + solver = get_model() + problem = first(solver.problems) +# problem.properties.finite_strain = true + problem.properties.geometric_stiffness = true + solver.properties.geometric_stiffness = true + solver(; debug=true) + @test isapprox(solver.properties.eigvals, [5/3, 2/3]) end @testset "test poisson problem modal analysis without tie" begin @@ -58,10 +63,6 @@ end 6 => [1.0, 3.0], 7 => [1.0, 9.0], 8 => [0.0, 9.0]) - T = Dict{Int64, Float64}() - for i=1:8 - T[i] = 0.0 - end el1 = Element(Quad4, [1, 2, 3, 4]) el2 = Element(Quad4, [4, 3, 7, 8]) el3 = Element(Seg2, [1, 2]) @@ -69,15 +70,12 @@ end update!([el1, el2, el3, el4], "geometry", X) update!([el1, el2], "density", 6.0) update!([el1, el2], "temperature thermal conductivity", 36.0) - #update!([el1, el2], "temperature", 0.0 => T) - update!([el1, el2], "temperature", T) update!([el3, el4], "temperature 1", 0.0) p1 = Problem(Heat, "combined body", 1) p1.properties.formulation = "2D" p2 = Problem(Dirichlet, "fixed ends", 1, "temperature") push!(p1, el1, el2) push!(p2, el3, el4) - solver = Solver(Modal) push!(solver, p1, p2) solver() @@ -94,10 +92,6 @@ end 6 => [1.0, 3.0], 7 => [1.0, 9.0], 8 => [0.0, 9.0]) - T = Dict{Int64, Float64}() - for i=1:8 - T[i] = 0.0 - end el1 = Element(Quad4, [1, 2, 3, 4]) el2 = Element(Quad4, [5, 6, 7, 8]) el3 = Element(Seg2, [1, 2]) @@ -105,8 +99,6 @@ end el5 = Element(Seg2, [3, 4]) el6 = Element(Seg2, [5, 6]) update!([el1, el2, el3, el4, el5, el6], "geometry", X) - #update!([el1, el2], "temperature", 0.0 => T) - update!([el1, el2], "temperature", T) update!([el1, el2], "density", 6.0) update!([el1, el2], "temperature thermal conductivity", 36.0) update!([el3, el4], "temperature 1", 0.0) diff --git a/test/test_mortar_2d_contact.jl b/test/test_mortar_2d_contact.jl index f09a2d5..b58867b 100644 --- a/test/test_mortar_2d_contact.jl +++ b/test/test_mortar_2d_contact.jl @@ -81,5 +81,6 @@ end la = slave("reaction force", [0.0], 0.0) info("u = $u, la = $la") @test isapprox(u, [-0.2, -0.15]) - @test isapprox(la, [0.0, 30.375]) + @test isapprox(la, [0.0, -30.375]) + # FIXME end diff --git a/test/test_mortar_2d_mesh_tie.jl b/test/test_mortar_2d_mesh_tie.jl index f89d5c4..fc38a8c 100644 --- a/test/test_mortar_2d_mesh_tie.jl +++ b/test/test_mortar_2d_mesh_tie.jl @@ -187,8 +187,9 @@ end slave_elements = get_slave_elements(interface) node_ids, la = get_nodal_vector(slave_elements, "reaction force", 0.0) for lai in la - @test isapprox(lai, [0.0, 10.0]) + @test isapprox(lai, [0.0, -10.0]) end + # FIXME end function JuliaFEM.get_mesh(::Type{Val{Symbol("curved 2d block splitted to upper and lower")}}) diff --git a/test/test_postprocess_xdmf.jl b/test/test_postprocess_xdmf.jl index edbfe85..05d6c5e 100644 --- a/test/test_postprocess_xdmf.jl +++ b/test/test_postprocess_xdmf.jl @@ -72,58 +72,6 @@ testdata = """\ """ -function test_write_to_xml() - nodes = Vector{Float64}[ - [0.0, 0.0, 0.0], - [1.0, 0.0, 0.0], - [0.0, 1.0, 0.0], - [0.0, 0.0, 1.0], - [0.5, 0.0, 0.0], - [0.5, 0.5, 0.0], - [0.0, 0.5, 0.0], - [0.0, 0.0, 0.5], - [0.5, 0.0, 0.5], - [0.0, 0.5, 0.5], - [1.0, 1.0, 1.0], - [2.0, 1.0, 1.0], - [1.0, 2.0, 1.0], - [1.0, 1.0, 2.0], - [1.5, 1.0, 1.0], - [1.5, 1.5, 1.0], - [1.0, 1.5, 1.0], - [1.0, 1.0, 1.5], - [1.5, 1.0, 1.5], - [1.0, 1.5, 1.5]] - - elements = [ - (:Tet10, [ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]) - (:Tet10, [11, 12, 13, 14, 15, 16, 17, 18, 19, 20])] - - displacement_field = nodes # same structure - xdoc, model = JuliaFEM.Postprocess.xdmf_new_model() - temporal_collection = JuliaFEM.Postprocess.xdmf_new_temporal_collection(model) - grid = JuliaFEM.Postprocess.xdmf_new_grid(temporal_collection; time=1) - JuliaFEM.Postprocess.xdmf_new_mesh!(grid, nodes, elements) - JuliaFEM.Postprocess.xdmf_new_nodal_field!(grid, "Displacement", displacement_field) - JuliaFEM.Postprocess.xdmf_save_model(xdoc, "/tmp/foo.xmf") - #info("exported data model: \n$(string(xdoc))") - #@test string(xdoc) == testdata - d1 = split(string(xdoc), "\n") -# d2 = split(testdata, "\n") - d2 = open(readlines, Pkg.dir("JuliaFEM")*"/test/testdata/quad_two_tet10.xmf") - println("comparing string") - for i in 1:length(d1) - println("d1: $(d1[i])") - println("d2: $(d2[i])") - #status = d1 == d2 ? "MATCHES" : "NO MATCH" - #info("line: $(d1[i]) $status") - #if d1 != d2 - # info("should be:\n$(d2[i])") - #end - d1 == d2 || error("No match") - end -end - @testset "write simple xmf file" begin X = Dict{Int64, Vector{Float64}}( 1 => [0.0, 0.0], diff --git a/test/test_problem.jl b/test/test_problem.jl index a93cf85..035d425 100644 --- a/test/test_problem.jl +++ b/test/test_problem.jl @@ -4,7 +4,7 @@ using JuliaFEM using JuliaFEM.Testing -@testset "test initialize field problem" begin +@testset "test initialize scalar field problem" begin el = Element(Seg2, [1, 2]) pr = Problem(Heat, 1) push!(pr, el) @@ -19,6 +19,20 @@ using JuliaFEM.Testing @test length(last(el, "temperature").data) == 2 end +@testset "test initialize vector field problem" begin + el = Element(Seg2, [1, 2]) + pr = Problem(Elasticity, 2) + push!(pr, el) + initialize!(pr) + @test haskey(el, "displacement") + @test length(el["displacement"]) == 1 + # this way we access to field at default time t=0.0, it's different than ^! + @test length(el("displacement")) == 2 + # length of single increment + @test length(el("displacement", 0.0)) == 2 + @test length(last(el, "displacement").data) == 2 +end + @testset "test initialize boundary problem" begin el = Element(Seg2, [1, 2]) pr = Problem(Dirichlet, "bc", 1, "temperature")