diff --git a/src/abaqus_reader.jl b/src/abaqus_reader.jl index 8c940c7..971b3f6 100644 --- a/src/abaqus_reader.jl +++ b/src/abaqus_reader.jl @@ -1,8 +1,6 @@ # This file is a part of JuliaFEM. # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md - - global handlers = Dict() """ @@ -43,7 +41,7 @@ end function parse_element_section(model, header, data) - Logging.debug("Parsing elements") + info("Parsing elements") eldims = Dict( "C3D10" => 10, "C3D4" => 4) @@ -58,14 +56,14 @@ function parse_element_section(model, header, data) elements = create_or_get(model, "elements") m = reshape(m, eldim+1, round(Int, length(m)/(eldim+1))) nel = size(m)[2] - Logging.debug("$nel elements found") + info("$nel elements found") for i=1:nel elements[m[1,i]] = m[2:end,i] end if "ELSET" in keys(header["options"]) elsets = create_or_get(model, "elsets") elset_name = header["options"]["ELSET"] - Logging.info("Creating ELSET $elset_name") + info("Creating ELSET $elset_name") elsets[elset_name] = Int64[] for i=1:nel push!(elsets[elset_name], m[1,i]) @@ -76,7 +74,7 @@ end function parse_nodeset_section(model, header, data) nset_name = header["options"]["NSET"] - Logging.debug("Creating node set $nset_name") + info("Creating node set $nset_name") m = matchall(r"[0-9]+", data) node_ids = map((s) -> parse(Int, s), m) nsets = create_or_get(model, "nsets") @@ -87,26 +85,25 @@ function parse_nodeset_section(model, header, data) end -function parse_abaqus(fid) +function parse_abaqus(fid::IOStream) model = Dict() - section = None - header = None - data = String[] - Logging.info("Registered handlers: $(keys(handlers))") + section = nothing + header = nothing + data = ASCIIString[] + info("Registered handlers: $(keys(handlers))") function process_section(section) - if section == None + if section == nothing return end if !(section in keys(handlers)) - Logging.info("Don't know what to do with data in section $section") - Logging.info("Skipping $(length(data)) bytes of unknown data") + info("Don't know what to do with data in section $section") + info("Skipping $(length(data)) bytes of unknown data") return end - Logging.debug("Starting to process") joined = join(data, "") handlers[section](model, header, strip(joined)) - data = [] + empty!(data) end line_idx = 0 @@ -115,16 +112,13 @@ function parse_abaqus(fid) continue end if startswith(line, "*") - Logging.debug("processing -- section") process_section(section) header = parse_header(line) - Logging.debug("Found ", header["section"], " section") section = header["section"] - Logging.debug("processing ++ section") continue end push!(data, line) - end + end process_section(section) return model end diff --git a/src/elasticity.jl b/src/elasticity.jl index af85b30..ed90847 100644 --- a/src/elasticity.jl +++ b/src/elasticity.jl @@ -94,6 +94,10 @@ function get_residual_vector(equation::ElasticityEquation, ip::IntegrationPoint, return vec(r) end + + + + ### Plane stress elasticity ### abstract PlaneElasticityProblem <: ElasticityProblem @@ -157,11 +161,79 @@ function get_residual_vector(equation::CPS2, ip::IntegrationPoint, time::Number; if haskey(element, "displacement traction force") T = basis("displacement traction force", ip, time) -# info("traction force = $T") -# info("basis = $(basis(ip, time))") r -= T*basis(ip, time) end return vec(r) end + +#= + +### 3d continuum elasticity ### + +type ContinuumElasticityProblem <: ElasticityProblem + unknown_field_name :: ASCIIString + unknown_field_dimension :: Int + equations :: Vector{ElasticityEquation} +end + +function ContinuumElasticityProblem(equations=[]) + return PlaneStressElasticityProblem("displacement", 3, equations) +end + +### Equations ### + +""" 4-node plane stress element. """ +type C3D10 <: ContinuumElasticityEquation + element :: Quad4 + integration_points :: Vector{IntegrationPoint} +end + +function Base.size(equation::CPS4) + return (2, 4) +end + +function Base.convert(::Type{PlaneStressElasticityEquation}, element::Quad4) + integration_points = get_integration_points(element) + if !haskey(element, "displacement") + element["displacement"] = 0.0 => [zeros(2) for i=1:4] + end + CPS4(element, integration_points) +end + +""" Boundary element for plane stress problem for surface loads. """ +type CPS2 <: PlaneStressElasticityEquation + element :: Seg2 + integration_points :: Vector{IntegrationPoint} +end + +function Base.size(equation::CPS2) + return (2, 2) +end + +function Base.convert(::Type{PlaneStressElasticityEquation}, element::Seg2) + integration_points = get_integration_points(element) + if !haskey(element, "displacement") + element["displacement"] = 0.0 => [zeros(2) for i=1:2] + end + CPS2(element, integration_points) +end + +function get_residual_vector(equation::CPS2, ip::IntegrationPoint, time::Number; variation=nothing) + + element = get_element(equation) + basis = get_basis(element) + + u = basis("displacement", ip, time, variation) + r = zeros(size(equation)) + + if haskey(element, "displacement traction force") + T = basis("displacement traction force", ip, time) + r -= T*basis(ip, time) + end + + return vec(r) +end + +=# diff --git a/src/integrate.jl b/src/integrate.jl index bd2cf24..97fed45 100644 --- a/src/integrate.jl +++ b/src/integrate.jl @@ -62,3 +62,20 @@ function get_integration_points(element::Seg3) return get_integration_points(element, Val{2}) end +### 3D elements + +function get_integration_points(element::Tet10, ::Type{Val{4}}) + a = .585410196624969 + b = .138196601125011 + w = .041666666666667 + integration_points = [ + IntegrationPoint([a, b, b], w), + IntegrationPoint([b, a, b], w), + IntegrationPoint([b, b, a], w), + IntegrationPoint([b, b, b], w)] +end + +function get_integration_points(element::Tet10) + return get_integration_points(element, Val{4}) +end +