diff --git a/Project.toml b/Project.toml index 0db850f..399039f 100644 --- a/Project.toml +++ b/Project.toml @@ -19,6 +19,7 @@ HeatTransfer = "4030f512-cedb-5907-ac7f-4ab05ad75ee7" InterfaceMechanics = "d2649932-d9e0-11e8-27ac-0193def238a7" LightXML = "9c8b4983-aa76-5018-a973-4c85ecc9e179" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +Logging = "56ddb016-857b-54e1-b83d-db4d58db5568" MortarContact2D = "048d6160-1a0b-53cd-a5b3-316946cc8d80" MortarContact2DAD = "c1673bdb-6aff-560b-99da-c78ea6da9af3" Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a" diff --git a/src/JuliaFEM.jl b/src/JuliaFEM.jl index 46cbad8..0a36479 100644 --- a/src/JuliaFEM.jl +++ b/src/JuliaFEM.jl @@ -113,6 +113,7 @@ import Base: getindex, setindex!, convert, length, size, isapprox, using SparseArrays, LinearAlgebra, Statistics using Reexport, ForwardDiff, LightXML, HDF5, Parameters +using Logging # For mesh readers using Tensors # For basis functions (Vec type) import Calculus # For symbolic differentiation in basis generation @@ -163,12 +164,15 @@ include("solvers/solvers_base.jl") # Base solver types include("analysis.jl") # Analysis and AbstractResultsWriter include("deprecated_fembase.jl") # Deprecated/legacy methods from FEMBase (length, size, etc.) +# Mesh readers (consolidated from AbaqusReader.jl and AsterReader.jl) +include("readers.jl") + using TimerOutputs export @timeit, print_timer # TODO: Consolidate these vendor packages later -# using AbaqusReader -# using AsterReader +# using AbaqusReader # Consolidated into src/readers.jl +# using AsterReader # Consolidated into src/readers.jl # @reexport using HeatTransfer include("problems_elasticity.jl") diff --git a/src/readers.jl b/src/readers.jl new file mode 100644 index 0000000..7ce0c0e --- /dev/null +++ b/src/readers.jl @@ -0,0 +1,15 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE +# +# Mesh readers consolidated from AbaqusReader.jl and AsterReader.jl + +# AbaqusReader - ABAQUS .inp file format +include("readers/keyword_register.jl") +include("readers/parse_mesh.jl") +include("readers/parse_model.jl") +include("readers/create_surface_elements.jl") +include("readers/abaqus_download.jl") + +# AsterReader - Code Aster .med file format +include("readers/read_aster_mesh.jl") +include("readers/read_aster_results.jl") diff --git a/src/readers/AbaqusReader.jl b/src/readers/AbaqusReader.jl new file mode 100644 index 0000000..bb54aaf --- /dev/null +++ b/src/readers/AbaqusReader.jl @@ -0,0 +1,17 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/AbaqusReader.jl/blob/master/LICENSE + +module AbaqusReader + +using Nullables + +include("parse_mesh.jl") +include("keyword_register.jl") +include("parse_model.jl") +include("create_surface_elements.jl") +include("abaqus_download.jl") + +export abaqus_read_mesh, abaqus_read_model, create_surface_elements +export abaqus_download + +end diff --git a/src/readers/AsterReader.jl b/src/readers/AsterReader.jl new file mode 100644 index 0000000..c22e02c --- /dev/null +++ b/src/readers/AsterReader.jl @@ -0,0 +1,25 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/AsterReader.jl/blob/master/LICENSE + +module AsterReader + +using HDF5 + +""" + parse_node_id(node_name) + +Return id number from node name. Usually the node name contains the id number, +i.e. N123 => 123 and so on. + +""" +function parse_node_id(node_name) + m = match(r"\d+", node_name) + m != nothing || error("Unable to parse id from node name $node_name.") + return tryparse(Int, m.match) +end + +include("read_aster_mesh.jl") +include("read_aster_results.jl") +export aster_read_mesh + +end diff --git a/src/readers/abaqus_download.jl b/src/readers/abaqus_download.jl new file mode 100644 index 0000000..5c2c4b1 --- /dev/null +++ b/src/readers/abaqus_download.jl @@ -0,0 +1,43 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/AbaqusReader.jl/blob/master/LICENSE + +""" + abaqus_download(model_name; dryrun=false) + +Download ABAQUS model from Internet. `model_name` is the name of the input +file. + +Given some model name from documentation, e.g., `et22sfse`, download that +file to local file system. This function uses environment variables to +determine the download url and place of storage. + +In order to use this function, one must set environment variable +`ABAQUS_DOWNLOAD_URL`, which determines a location where to download. For +example, if the path to model is `https://domain.com/v6.14/books/eif/et22sfse.inp`, +`ABAQUS_DOWNLOAD_URL` will be the basename of that path, i.e., +`https://domain.com/v6.14/books/eif`. + +By default, the model will be downloaded to current directory. If that is not +desired, one can set another environment variable `ABAQUS_DOWNLOAD_DIR`, and +in that case the file will be downloaded to that directory. + +Function call will return full path to downloaded file or nothing, if download +is failing because of missing environment variable `ABAQUS_DOWNLOAD_DIR`. +""" +function abaqus_download(model_name, env=ENV; dryrun=false) + path = get(env, "ABAQUS_DOWNLOAD_DIR", "") + fn = joinpath(path, model_name) + if isfile(fn) # already downloaded + return fn + end + if !haskey(env, "ABAQUS_DOWNLOAD_URL") + error("ABAQUS input file $fn not found and `ABAQUS_DOWNLOAD_URL` not ", + "set, unable to download file. To enable automatic model ", + "downloading, set url to models to environment variable + `ABAQUS_DOWNLOAD_URL`") + end + url = joinpath(env["ABAQUS_DOWNLOAD_URL"], model_name) + @debug("Downloading model $model_name to $fn") + dryrun || download(url, fn) + return fn +end diff --git a/src/readers/create_surface_elements.jl b/src/readers/create_surface_elements.jl new file mode 100644 index 0000000..b2a4260 --- /dev/null +++ b/src/readers/create_surface_elements.jl @@ -0,0 +1,68 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/AbaqusReader.jl/blob/master/LICENSE + +""" + element_mapping + +This mapping table contains information what node ids locally match +each side of element. +""" +const element_mapping = Dict( + :Tet4 => Dict( + :S1 => (:Tri3, [1, 3, 2]), + :S2 => (:Tri3, [1, 2, 4]), + :S3 => (:Tri3, [2, 3, 4]), + :S4 => (:Tri3, [1, 4, 3])), + :Tet10 => Dict( + :S1 => (:Tri6, [1, 3, 2, 7, 6, 5]), + :S2 => (:Tri6, [1, 2, 4, 5, 9, 8]), + :S3 => (:Tri6, [2, 3, 4, 6, 10, 9]), + :S4 => (:Tri6, [1, 4, 3, 8, 10, 7])), + :Hex8 => Dict( + :S1 => (:Quad4, [1, 2, 3, 4]), + :S2 => (:Quad4, [5, 8, 7, 6]), + :S3 => (:Quad4, [1, 5, 6, 2]), + :S4 => (:Quad4, [2, 6, 7, 3]), + :S5 => (:Quad4, [3, 7, 8, 4]), + :S6 => (:Quad4, [4, 8, 5, 1]))) + +""" Given element code, element side and global connectivity, determine boundary +element. E.g. for Tet4 we have 4 sides S1..S4 and boundary element is of type Tri3. +""" +function create_surface_element(element_type::Symbol, element_side::Symbol, + element_connectivity::Vector{Int}) + + if !haskey(element_mapping, element_type) + error("Unable to find surface element for element of type ", + "$element_type for side $element_side, update element ", + "mapping table.") + 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, update ", + "element mapping table.") + end + + surfel, surfel_lconn = element_mapping[element_type][element_side] + surfel_gconn = element_connectivity[surfel_lconn] + return surfel, surfel_gconn +end + +""" + create_surface_elements(mesh, surface_name) + +Create surface elements for `surface` using mesh `mesh`. +Mesh can be obtained by using `abaqus_read_mesh`. +""" +function create_surface_elements(mesh::Dict, surface_name::String) + surface = mesh["surface_sets"][surface_name] + elements = mesh["elements"] + eltypes = mesh["element_types"] + result = Tuple{Symbol, Vector{Int}}[] + for (elid, side) in surface + surface_element = create_surface_element(eltypes[elid], side, elements[elid]) + push!(result, surface_element) + end + return result +end diff --git a/src/readers/keyword_register.jl b/src/readers/keyword_register.jl new file mode 100644 index 0000000..62b3755 --- /dev/null +++ b/src/readers/keyword_register.jl @@ -0,0 +1,24 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/AbaqusReader.jl/blob/master/LICENSE + +global __register__ = Set{String}() + +""" + register_abaqus_keyword(keyword::String) + +Add ABAQUS keyword `s` to register. That is, after registration every time +keyword show up in `.inp` file a new section is started +""" +function register_abaqus_keyword(keyword::String) + push!(__register__, keyword) + return Type{Val{Symbol(keyword)}} +end + +""" + is_abaqus_keyword_registered(keyword::String) + +Return true/false is ABAQUS keyword registered. +""" +function is_abaqus_keyword_registered(keyword::String) + return keyword in __register__ +end diff --git a/src/readers/parse_mesh.jl b/src/readers/parse_mesh.jl new file mode 100644 index 0000000..90bab5d --- /dev/null +++ b/src/readers/parse_mesh.jl @@ -0,0 +1,305 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/AbaqusReader.jl/blob/master/LICENSE + +import Base.parse +# using Logging # Already imported in JuliaFEM.jl + +# Define element type and number of nodes in element +element_has_nodes(::Type{Val{:C3D4}}) = 4 +element_has_type(::Type{Val{:C3D4}}) = :Tet4 + +element_has_nodes(::Type{Val{:C3D6}}) = 6 +element_has_type(::Type{Val{:C3D6}}) = :Wedge6 + +element_has_nodes(::Type{Val{:C3D4H}}) = 4 +element_has_type(::Type{Val{:C3D4H}}) = :Tet4 + +element_has_nodes(::Type{Val{:C3D8}}) = 8 +element_has_type(::Type{Val{:C3D8}}) = :Hex8 + +element_has_nodes(::Type{Val{:C3D8R}}) = 8 +element_has_type(::Type{Val{:C3D8R}}) = :Hex8 + +element_has_nodes(::Type{Val{:COH3D8}}) = 8 +element_has_type(::Type{Val{:COH3D8}}) = :Hex8 + +element_has_nodes(::Type{Val{:C3D10}}) = 10 +element_has_type(::Type{Val{:C3D10}}) = :Tet10 + +element_has_nodes(::Type{Val{:C3D10H}}) = 10 +element_has_type(::Type{Val{:C3D10H}}) = :Tet10 + +element_has_nodes(::Type{Val{:C3D20}}) = 20 +element_has_type(::Type{Val{:C3D20}}) = :Hex20 + +element_has_nodes(::Type{Val{:C3D20E}}) = 20 + +element_has_nodes(::Type{Val{:S3}}) = 3 +element_has_type(::Type{Val{:S3}}) = :Tri3 + +element_has_nodes(::Type{Val{:CPS3}}) = 3 +element_has_type(::Type{Val{:CPS3}}) = :CPS3 + +element_has_nodes(::Type{Val{:STRI65}}) = 6 +element_has_type(::Type{Val{:STRI65}}) = :Tri6 + +element_has_nodes(::Type{Val{:CPS4}}) = 4 +element_has_type(::Type{Val{:CPS4}}) = :Quad4 + +element_has_nodes(::Type{Val{:CPS4R}}) = 4 +element_has_type(::Type{Val{:CPS4R}}) = :Quad4 + +element_has_nodes(::Type{Val{:T2D2}}) = 2 +element_has_type(::Type{Val{:T2D2}}) = :Seg2 + +element_has_nodes(::Type{Val{:T3D2}}) = 2 +element_has_type(::Type{Val{:T3D2}}) = :Seg2 + +element_has_nodes(::Type{Val{:B33}}) = 2 +element_has_type(::Type{Val{:B33}}) = :Seg2 + +"""Checks for a comment or empty line + +Function return true, if line starts with comment character "**" +or has length of 0 +""" +function empty_or_comment_line(line::T) where {T<:AbstractString} + startswith(line, "**") || (length(line) == 0) +end + +"""Match words from both sides of '=' character +""" +function matchset(definition) + regexp = r"([\w\_\-]+[ ]*=[ ]*[\w\_\-]+)" + collect(m.match for m = eachmatch(regexp, definition)) +end + +"""Parse string to get set type and name +""" +function parse_definition(definition) + set_defs = Dict() + set_definition = matchset(definition) + set_definition == nothing && return nothing + for x in set_definition + name, vals = map(strip, split(x, "=")) + set_defs[lowercase(name)] = vals + end + set_defs +end + +"""Parse all the numbers from string +""" +function parse_numbers(line, type_::Type{T})::Vector{T} where {T} + regexp = r"[0-9]+" + matches = collect((m.match for m = eachmatch(regexp, line))) + map(x -> Base.parse(type_, x), matches) +end + +"""Add set to model, if set exists +""" +function add_set!(model, definition, model_key, abaqus_key, ids) + has_set_def = parse_definition(definition) + if haskey(has_set_def, "elset") + set_name = has_set_def[abaqus_key] + @debug("Adding $abaqus_key: $set_name") + model[model_key][set_name] = ids + end +end + +"""Parse nodes from the lines +""" +function parse_section(model, lines, ::Symbol, idx_start, idx_end, ::Type{Val{:NODE}}) + nnodes = 0 + ids = Int[] + definition = lines[idx_start] + for line in lines[idx_start+1:idx_end] + if !(empty_or_comment_line(line)) + m = collect((m.match for m = eachmatch(r"[-0-9.eE+]+", line))) + node_id = parse(Int, m[1]) + coords = parse.(Float64, m[2:end]) + model["nodes"][node_id] = coords + push!(ids, node_id) + nnodes += 1 + end + end + @debug("$nnodes nodes found") + add_set!(model, definition, "node_sets", "nset", ids) +end + +"""Custon regex to find match from string. Index used if there are multiple matches +""" +function regex_match(regex_str, line, idx) + return match(regex_str, line).captures[idx] +end + +"""Custom list iterator + +Simple iterator for comsuming element list. Depending +on the used element, connectivity nodes might be listed +in multiple lines, which is why iterator is used to handle +this problem. +""" +function consumeList(arr, start, stop) + idx = start - 1 + function _it() + idx += 1 + if idx > stop + return nothing + end + arr[idx] + end + _it +end + +"""Parse elements from input lines + +Reads element ids and their connectivity nodes from input lines. +If elset definition exists, also adds the set to model. +""" +function parse_section(model, lines, ::Symbol, idx_start, idx_end, ::Type{Val{:ELEMENT}}) + ids = Int[] + definition = lines[idx_start] + regexp = r"TYPE=([\w\-\_]+)"i + m = match(regexp, definition) + m == nothing && error("Could not match regexp $regexp to line $definition") + element_type = uppercase(m[1]) + eltype_sym = Symbol(element_type) + eltype_nodes = element_has_nodes(Val{eltype_sym}) + element_type = element_has_type(Val{eltype_sym}) + @debug("Parsing elements. Type: $(m[1]). Topology: $(element_type)") + list_iterator = consumeList(lines, idx_start + 1, idx_end) + line = list_iterator() + while line != nothing + numbers = parse_numbers(line, Int) + if !(empty_or_comment_line(line)) + id = numbers[1] + push!(ids, id) + connectivity = numbers[2:end] + while length(connectivity) != eltype_nodes + @assert length(connectivity) < eltype_nodes + line = list_iterator() + numbers = parse_numbers(line, Int) + push!(connectivity, numbers...) + end + model["elements"][id] = connectivity + model["element_types"][id] = element_type + end + line = list_iterator() + end + add_set!(model, definition, "element_sets", "elset", ids) +end + +"""Parse node and elementset from input lines +""" +function parse_section(model, lines, key, idx_start, idx_end, ::Union{Type{Val{:NSET}}, + Type{Val{:ELSET}}}) + data = Int[] + set_regex_string = Dict(:NSET => r"((?<=NSET=)([\w\-\_]+)|(?<=NSET=\")([\w\-\_\ ]+)(?=\"))"i, + :ELSET => r"((?<=ELSET=)([\w\-\_]+)|(?<=ELSET=\")([\w\-\_\ ]+)(?=\"))"i) + selected_set = key == :NSET ? "node_sets" : "element_sets" + definition = lines[idx_start] + regex_string = set_regex_string[key] + set_name = regex_match(regex_string, definition, 1) + @debug("Creating $(lowercase(string(key))) $set_name") + + if endswith(strip(uppercase(definition)), "GENERATE") + line = lines[idx_start+1] + first_id, last_id, step_ = parse_numbers(line, Int) + set_ids = collect(first_id:step_:last_id) + push!(data, set_ids...) + else + for line in lines[idx_start+1:idx_end] + if !(empty_or_comment_line(line)) + set_ids = parse_numbers(line, Int)::Vector{Int} + push!(data, set_ids...) + end + end + end + model[selected_set][set_name] = data +end + +"""Parse SURFACE keyword +""" +function parse_section(model, lines, ::Symbol, idx_start, idx_end, ::Type{Val{:SURFACE}}) + data = Vector{Tuple{Int,Symbol}}() + definition = lines[idx_start] + + has_set_def = parse_definition(definition) + has_set_def != nothing || return + set_type = get(has_set_def, "type", "UNKNOWN") + set_name = has_set_def["name"] + + for line in lines[idx_start+1:idx_end] + empty_or_comment_line(line) && continue + m = match(r"(?P\d+),.*(?PS\d+).*", line) + element_id = parse(Int, m[:element_id]) + element_side = Symbol(m[:element_side]) + push!(data, (element_id, element_side)) + end + model["surface_types"][set_name] = Symbol(set_type) + model["surface_sets"][set_name] = data + return +end + +"""Find lines, which contain keywords, for example "*NODE" +""" +function find_keywords(lines) + indexes = Int[] + for (idx, line) in enumerate(lines) + if startswith(line, "*") && !startswith(line, "**") + push!(indexes, idx) + end + end + return indexes +end + +"""Main function for parsing Abaqus input file. + +Function parses Abaqus input file and generates a dictionary of +all the available keywords. +""" +function parse_abaqus(fid::IOStream, verbose::Bool) + model = Dict{String,Dict}() + model["nodes"] = Dict{Int,Vector{Float64}}() + model["node_sets"] = Dict{String,Vector{Int}}() + model["elements"] = Dict{Int,Vector{Int}}() + model["element_types"] = Dict{Int,Symbol}() + model["element_sets"] = Dict{String,Vector{Int}}() + model["surface_sets"] = Dict{String,Vector{Tuple{Int,Symbol}}}() + model["surface_types"] = Dict{String,Symbol}() + keyword_sym::Symbol = :none + + lines = readlines(fid) + keyword_indexes = find_keywords(lines) + nkeyword_indexes = length(keyword_indexes) + push!(keyword_indexes, length(lines) + 1) + idx_start = keyword_indexes[1] + + for idx_end in keyword_indexes[2:end] + keyword_line = strip(uppercase(lines[idx_start])) + keyword = strip(regex_match(r"\s*([\w ]+)", keyword_line, 1)) + k_sym = Symbol(keyword) + args = Tuple{Dict,Vector{Int},Symbol,Int,Int,Type{Val{k_sym}}} + if hasmethod(parse_section, args) + parse_section(model, lines, k_sym, idx_start, idx_end - 1, Val{k_sym}) + else + if verbose + @warn("Unknown section: '$(keyword)'") + end + end + idx_start = idx_end + end + return model +end + +""" + abaqus_read_mesh(fn::String) + +Read ABAQUS mesh from file `fn`. Returns a dict with elements, nodes, +element sets, node sets and other topologically imporant things, but +not the actual model with boundary conditions, load steps and so on. +""" +function abaqus_read_mesh(fn::String; kwargs...) + verbose = get(kwargs, :verbose, true) + return parse_abaqus(open(fn), verbose) +end diff --git a/src/readers/parse_model.jl b/src/readers/parse_model.jl new file mode 100644 index 0000000..825f97c --- /dev/null +++ b/src/readers/parse_model.jl @@ -0,0 +1,347 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/AbaqusReader.jl/blob/master/LICENSE + +import Base.parse +import Base: getindex, length + +### Model definitions for ABAQUS data model + +abstract type AbstractMaterial end +abstract type AbstractMaterialProperty end +abstract type AbstractProperty end +abstract type AbstractStep end +abstract type AbstractBoundaryCondition end +abstract type AbstractOutputRequest end + +mutable struct Mesh + nodes :: Dict{Int, Vector{Float64}} + node_sets :: Dict{String, Vector{Int}} + elements :: Dict{Int, Vector{Int}} + element_types :: Dict{Int, Symbol} + element_sets :: Dict{String, Vector{Int}} + surface_sets :: Dict{String, Vector{Tuple{Int, Symbol}}} + surface_types :: Dict{String, Symbol} +end + +function Mesh(d::Dict{String, Dict}) + return Mesh(d["nodes"], d["node_sets"], d["elements"], + d["element_types"], d["element_sets"], + d["surface_sets"], d["surface_types"]) +end + +mutable struct Model + path :: String + name :: String + mesh :: Mesh + materials :: Dict{Symbol, AbstractMaterial} + properties :: Vector{AbstractProperty} + boundary_conditions :: Vector{AbstractBoundaryCondition} + steps :: Vector{AbstractStep} +end + +mutable struct SolidSection <: AbstractProperty + element_set :: Symbol + material_name :: Symbol +end + +mutable struct Material <: AbstractMaterial + name :: Symbol + properties :: Vector{AbstractMaterialProperty} +end + +mutable struct Elastic <: AbstractMaterialProperty + E :: Float64 + nu :: Float64 +end + +mutable struct Step <: AbstractStep + kind :: Union{Symbol, Nothing} # STATIC, ... (was Nullable{Symbol} in Julia 0.x) + boundary_conditions :: Vector{AbstractBoundaryCondition} + output_requests :: Vector{AbstractOutputRequest} +end + +mutable struct BoundaryCondition <: AbstractBoundaryCondition + kind :: Symbol # BOUNDARY, CLOAD, DLOAD, DSLOAD, ... + data :: Vector + options :: Dict +end + +mutable struct OutputRequest <: AbstractOutputRequest + kind :: Symbol # NODE, EL, SECTION, ... + data :: Vector + options :: Dict + target :: Symbol # PRINT, FILE +end + +### Utility functions to parse ABAQUS .inp file to data model + +mutable struct Keyword + name :: String + options :: Vector{Union{String, Pair}} +end + +mutable struct AbaqusReaderState + section :: Union{Keyword, Nothing} # was Nullable{Keyword} + material :: Union{AbstractMaterial, Nothing} # was Nullable + property :: Union{AbstractProperty, Nothing} # was Nullable + step :: Union{AbstractStep, Nothing} # was Nullable + data :: Vector{String} +end + +function get_data(state::AbaqusReaderState) + data = [] + for row in state.data + row = strip(row, [' ', ',']) + col = split(row, ',') + col = map(Meta.parse, col) + push!(data, col) + end + return data +end + +function get_options(state::AbaqusReaderState) + return Dict(state.section.options) +end + +function get_option(state::AbaqusReaderState, what::String) + return get_options(state)[what] +end + +function length(state::AbaqusReaderState) + return length(state.data) +end + +function is_comment(line) + return startswith(line, "**") +end + +function is_keyword(line) + return startswith(line, "*") && !is_comment(line) +end + +function parse_keyword(line; uppercase_keyword=true) + args = split(line, ",") + args = map(String, map(strip, args)) + keyword_name = strip(args[1], '*') + if uppercase_keyword + keyword_name = uppercase(keyword_name) + end + keyword = Keyword(keyword_name, []) + for option in args[2:end] + pair = map(String, split(option, "=")) + if uppercase_keyword + pair[1] = uppercase(pair[1]) + end + if length(pair) == 1 + push!(keyword.options, pair[1]) + elseif length(pair) == 2 + push!(keyword.options, pair[1] => pair[2]) + else + error("Keyword failure: $line, $option, $pair") + end + end + return keyword +end + +function is_new_section(line) + is_keyword(line) || return false + section = parse_keyword(line) + is_abaqus_keyword_registered(section.name) || return false + return true +end + +# open_section! is called right after keyword is found +function open_section! end + +# close_section! is called at the end or section or before new keyword +function close_section! end + +function maybe_close_section!(model, state) + global close_section! + isnull(state.section) && return + section_name = state.section.name + @debug("Close section: $section_name") + args = Tuple{Model, AbaqusReaderState, Type{Val{Symbol(section_name)}}} + if hasmethod(close_section!, args) + close_section!(model, state, Val{Symbol(section_name)}) + else + @debug("no close_section! found for $section_name") + end + state.section = nothing +end + +function maybe_open_section!(model, state) + global open_section! + section_name = state.section.name + section_options = state.section.options + @debug("New section: $section_name with options $section_options") + args = Tuple{Model, AbaqusReaderState, Type{Val{Symbol(section_name)}}} + if hasmethod(open_section!, args) + open_section!(model, state, Val{Symbol(section_name)}) + else + @warn("no open_section! found for $section_name") + end +end + +function new_section!(model, state, line::String) + maybe_close_section!(model, state) + state.data = [] + state.section = parse_keyword(line) + maybe_open_section!(model, state) +end + +function process_line!(model, state, line::String) + if isnull(state.section) + @debug("section = nothing! line = $line") + return + end + if is_keyword(line) + @warn("missing keyword? line = $line") + # close section, this is probably keyword and collecting data should stop. + maybe_close_section!(model, state) + return + end + push!(state.data, line) +end + +""" + abaqus_read_model(filename::String) + +Read ABAQUS model from file. Include also boundary conditions, load steps +and so on. If only mesh is needed, it's better to use `abaqus_read_mesh` +insted. +""" +function abaqus_read_model(fn::String) + + model_path = dirname(fn) + model_name = first(splitext(basename(fn))) + mesh = Mesh(open(parse_abaqus, fn)) + materials = Dict() + model = Model(model_path, model_name, mesh, materials, [], [], []) + + state = AbaqusReaderState(nothing, nothing, nothing, nothing, []) + + fid = open(fn) + for line in eachline(fid) + line = convert(String, strip(line)) + is_comment(line) && continue + if is_new_section(line) + new_section!(model, state, line) + else + process_line!(model, state, line) + end + end + close(fid) + maybe_close_section!(model, state) + + return model +end + +### Code to parse ABAQUS .inp to data model + +# add here only keywords when planning to define open_section! and/or +# close_section!, i.e. actually parse keyword to model. little bit of +# magic is happening here, but after calling macro there is typealias +# defined i.e. typealias SOLID_SECTION Type{Val{Symbol("SOLID_SECTION")}} +# and also is_keyword_registered("SOLID SECTION") returns true after +# registration, also notice underscoring + +SOLID_SECTION = register_abaqus_keyword("SOLID SECTION") + +MATERIAL = register_abaqus_keyword("MATERIAL") +ELASTIC = register_abaqus_keyword("ELASTIC") + +STEP = register_abaqus_keyword("STEP") +STATIC = register_abaqus_keyword("STATIC") +END_STEP = register_abaqus_keyword("END STEP") + +BOUNDARY = register_abaqus_keyword("BOUNDARY") +CLOAD = register_abaqus_keyword("CLOAD") +DLOAD = register_abaqus_keyword("DLOAD") +DSLOAD = register_abaqus_keyword("DSLOAD") +const BOUNDARY_CONDITIONS = Union{BOUNDARY, CLOAD, DLOAD, DSLOAD} + +NODE_PRINT = register_abaqus_keyword("NODE PRINT") +EL_PRINT = register_abaqus_keyword("EL PRINT") +SECTION_PRINT = register_abaqus_keyword("SECTION PRINT") +const OUTPUT_REQUESTS = Union{NODE_PRINT, EL_PRINT, SECTION_PRINT} + +## Properties + +function open_section!(model, state, ::SOLID_SECTION) + element_set = Symbol(get_option(state, "ELSET")) + material_name = Symbol(get_option(state, "MATERIAL")) + property = SolidSection(element_set, material_name) + state.property = property + push!(model.properties, property) +end + +function close_section!(model, state, ::SOLID_SECTION) + name = model.name + state.property = nothing +end + +## Materials + +function open_section!(model, state, ::MATERIAL) + material_name = Symbol(get_option(state, "NAME")) + material = Material(material_name, []) + state.material = material + model.materials[material_name] = material +end + +function close_section!(model, state, ::ELASTIC) + name = model.name + @assert length(state) == 1 + E, nu = first(get_data(state)) + material_property = Elastic(E, nu) + material = state.material + push!(material.properties, material_property) +end + +## Steps + +function open_section!(model, state, ::STEP) + name = model.name + step_ = Step(nothing, Vector(), Vector()) + state.step = step_ + push!(model.steps, step_) +end + +function open_section!(model, state, ::STATIC) + name = model.name + isnull(state.step) && error("*STATIC outside *STEP ?") + state.step.kind = :STATIC +end + +function open_section!(model, state, ::END_STEP) + name = model.name + state.step = nothing +end + +## Steps -- boundary conditions + +function close_section!(model, state, ::BOUNDARY_CONDITIONS) + kind = Symbol(state.section.name) + data = get_data(state) + options = get_options(state) + bc = BoundaryCondition(kind, data, options) + if isnull(state.step) + push!(model.boundary_conditions, bc) + else + step_ = state.step + push!(step_.boundary_conditions, bc) + end +end + +## Steps -- output requests + +function close_section!(model, state, ::OUTPUT_REQUESTS) + name = model.name + kind, target = map(Meta.parse, split(state.section.name, " ")) + data = get_data(state) + options = get_options(state) + request = OutputRequest(Symbol(kind), data, options, Symbol(target)) + step_ = state.step + push!(step_.output_requests, request) +end diff --git a/src/readers/read_aster_mesh.jl b/src/readers/read_aster_mesh.jl new file mode 100644 index 0000000..6909caf --- /dev/null +++ b/src/readers/read_aster_mesh.jl @@ -0,0 +1,222 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/AsterReader.jl/blob/master/LICENSE + +function aster_parse_nodes(section; strip_characters=true) + nodes = Dict{Any, Vector{Float64}}() + has_started = false + for line in split(section, '\n') + m = collect((string(m.match) for m in eachmatch(r"[\w.-]+", line))) + if length(m) == 1 + if (m[1] == "COOR_2D") || (m[1] == "COOR_3D") + has_started = true + continue + end + if m[1] == "FINSF" + break + end + end + if !has_started # we have not found COOR_2D or COOR_3D yet. + continue + end + to(T, x) = map(x -> parse(T, x), x) + if length(m) == 4 + if strip_characters + nodes[parse_node_id(m[1])] = to(Float64, m[2:end]) + else + nodes[nid] = to(Float64, m[2:end]) + end + end + end + return nodes +end + + +""" Code Aster binary file (.med). """ +mutable struct MEDFile + data :: Dict +end + +function MEDFile(fn::String) + return MEDFile(h5read(fn, "/")) +end + +function get_mesh_names(med::MEDFile) + return sort(collect(keys(med.data["FAS"]))) +end + +""" Convert vector of Int8 to ASCII string. """ +function to_ascii(data::Vector{Int8}) + return ascii(unsafe_string(pointer(convert(Vector{UInt8}, data)))) +end + +function get_mesh(med::MEDFile, mesh_name::String) + if !haskey(med.data["FAS"], mesh_name) + @warn("Mesh $mesh_name not found from med file.") + meshes = get_mesh_names(med) + all_meshes = join(meshes, ", ") + @warn("Available meshes: $all_meshes") + error("Mesh $mesh_name not found.") + end + return med.data["FAS"][mesh_name] +end + +""" Return node sets from med file. + +Notes +----- +One node set id can have multiple names. + +""" +function get_node_sets(med::MEDFile, mesh_name::String)::Dict{Int64, Vector{String}} + mesh = get_mesh(med, mesh_name) + node_sets = Dict{Int64, Vector{String}}(0 => ["OTHER"]) + if !haskey(mesh, "NOEUD") + return node_sets + end + for (k, v) in mesh["NOEUD"] + nset_id = parse(Int, split(k, "_")[2]) + node_sets[nset_id] = collect(to_ascii(d) for d in v["GRO"]["NOM"]) + end + return node_sets +end + +""" + get_element_sets(med, mesh_name) + +Return element sets from med file. Return type is a dictionary, where the key is +the element set id number (integer) and value is a vector of strings, containing +human-readable name for element set. + +# Notes + +One element set id can have multiple names. + +""" +function get_element_sets(med::MEDFile, mesh_name::String)::Dict{Int64, Vector{String}} + mesh = get_mesh(med, mesh_name) + element_sets = Dict{Int64, Vector{String}}() + if !haskey(mesh, "ELEME") + return element_sets + end + elset_keys = sort(collect(keys(mesh["ELEME"]))) + for (i, k) in enumerate(elset_keys) + v = mesh["ELEME"][k] + if startswith(k, "FAM") # exported from Salome + elset_id = parse(Int, split(k, '_')[2]) + else # exported from Gmsh + elset_id = -i + end + if !isempty(v) + element_sets[elset_id] = map(strip, collect(to_ascii(d) for d in v["GRO"]["NOM"])) + else + element_sets[elset_id] = [""] + end + end + return element_sets +end + +function get_nodes(med::MEDFile, nsets::Dict{Int, Vector{String}}, mesh_name::String) + increments = keys(med.data["ENS_MAA"][mesh_name]) + @assert length(increments) == 1 + increment = first(increments) + nodes = med.data["ENS_MAA"][mesh_name][increment]["NOE"] + nset_ids = nodes["FAM"] + nnodes = length(nset_ids) + node_ids = get(nodes, "NUM", collect(1:nnodes)) + node_coords = nodes["COO"] + dim = round(Int, length(node_coords)/nnodes) + node_coords = reshape(node_coords, nnodes, dim)' + d = Dict{Int64}{Tuple{Vector{String}, Vector{Float64}}}() + for i=1:nnodes + nset = nsets[nset_ids[i]] + d[node_ids[i]] = (nset, node_coords[:, i]) + end + return d +end + +function get_connectivity(med::MEDFile, elsets::Dict{Int64, Vector{String}}, mesh_name::String) + if !haskey(elsets, 0) + elsets[0] = ["OTHER"] + end + increments = keys(med.data["ENS_MAA"][mesh_name]) + @assert length(increments) == 1 + increment = first(increments) + all_elements = med.data["ENS_MAA"][mesh_name][increment]["MAI"] + d = Dict{Int64, Tuple{Symbol, Vector{String}, Vector{Int64}}}() + for eltype in keys(all_elements) + elements = all_elements[eltype] + elset_ids = elements["FAM"] + nelements = length(elset_ids) + element_ids = get(elements, "NUM", collect(1:nelements)) + element_connectivity = elements["NOD"] + element_dim = round(Int, length(element_connectivity)/nelements) + element_connectivity = reshape(element_connectivity, nelements, element_dim)' + for i=1:nelements + eltype = Symbol(eltype) + elco = element_connectivity[:, i] + elset = elsets[elset_ids[i]] + d[element_ids[i]] = (eltype, elset, elco) + end + end + return d +end + +""" + aster_read_mesh(filename, mesh_name=nothing) + +Parse code aster .med file and return mesh data in a dictionary. + +Dictionary contains additional dictionaries `nodes`, `node_sets`, `elements`, +`element_sets`, `element_types`, `surface_sets` and `surface_types`. + +If mesh file contains several meshes, one must provide the mesh name as +additional argument or expcetion will be thrown. + +""" +function aster_read_mesh(filename::String, mesh_name=nothing) + isfile(filename) || error("Cannot read mesh file $filename: file not found.") + aster_read_mesh_(MEDFile(filename), mesh_name) +end + +function aster_read_mesh_(med::MEDFile, mesh_name=nothing) + mesh_names = get_mesh_names(med) + all_meshes = join(mesh_names, ", ") + if mesh_name == nothing + length(mesh_names) == 1 || error("several meshes found from med, pick one: $all_meshes") + mesh_name = mesh_names[1] + else + mesh_name in mesh_names || error("Mesh $mesh_name not found from mesh file $fn. Available meshes: $all_meshes") + end + + mesh = Dict{String, Dict}() + mesh["nodes"] = Dict{Int, Vector{Float64}}() + mesh["node_sets"] = Dict{String, Vector{Int}}() + mesh["elements"] = Dict{Int, Vector{Int}}() + mesh["element_types"] = Dict{Int, Symbol}() + mesh["element_sets"] = Dict{String, Vector{Int}}() + mesh["surface_sets"] = Dict{String, Vector{Tuple{Int, Symbol}}}() + mesh["surface_types"] = Dict{String, Symbol}() + + elsets = get_element_sets(med, mesh_name) + nsets = get_node_sets(med, mesh_name) + for (nid, (nset_, coords)) in get_nodes(med, nsets, mesh_name) + mesh["nodes"][nid] = coords + for nset in nset_ + if !haskey(mesh["node_sets"], nset) + mesh["node_sets"][nset] = [] + end + push!(mesh["node_sets"][nset], nid) + end + end + for (elid, (eltyp, elset_, elcon)) in get_connectivity(med, elsets, mesh_name) + mesh["elements"][elid] = elcon + mesh["element_types"][elid] = eltyp + for elset in elset_ + if !haskey(mesh["element_sets"], elset) + mesh["element_sets"][elset] = [] + end + push!(mesh["element_sets"][elset], elid) + end + end + return mesh +end diff --git a/src/readers/read_aster_results.jl b/src/readers/read_aster_results.jl new file mode 100644 index 0000000..ba84f44 --- /dev/null +++ b/src/readers/read_aster_results.jl @@ -0,0 +1,60 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/AsterReader.jl/blob/master/LICENSE + +""" Code Aster result file (.rmed). """ +mutable struct RMEDFile + data :: Dict +end + +function RMEDFile(fn::String) + return RMEDFile(h5read(fn, "/")) +end + +""" Return nodes from result med file. """ +function aster_read_nodes(rmed::RMEDFile) + increments = keys(rmed.data["ENS_MAA"]["MAIL"]) + @assert length(increments) == 1 + increment = first(increments) + nodes = rmed.data["ENS_MAA"]["MAIL"][increment]["NOE"] + node_names = nodes["NOM"] + node_coords = nodes["COO"] + nnodes = length(node_names) + dim = round(Int, length(node_coords)/nnodes) + node_coords = reshape(node_coords, nnodes, dim)' + stripper(node_name) = strip(ascii(unsafe_string(pointer(convert(Vector{UInt8}, node_name))))) + node_names = map(stripper, node_names) + node_ids = map(parse_node_id, node_names) + nodes = Dict(j => node_coords[:,j] for j in node_ids) + return nodes +end + +""" Read nodal field from rmed file. """ +function aster_read_data(rmed::RMEDFile, field_name; field_type=:NODE, + info_fields=true, node_ids=nothing) + + if occursin("ELGA", field_name) + field_type = :GAUSS + end + + if node_ids == nothing + nodes = aster_read_nodes(rmed) + node_ids = sort(collect(keys(nodes))) + end + + if info_fields + field_names = keys(rmed.data["CHA"]) + all_fields = join(field_names, ", ") + @info("results: $all_fields") + end + + chdata = rmed.data["CHA"]["RESU____$field_name"] + @assert length(chdata) == 1 + increment = chdata[first(keys(chdata))] + if field_type == :NODE + data = increment["NOE"]["MED_NO_PROFILE_INTERNAL"]["CO"] + results = Dict(j => data[j] for j in node_ids) + else + error("Unable to read result of type $field_type: not implemented") + end + return results +end diff --git a/src/sparse/sparse.jl b/src/sparse/sparse.jl index cf0a812..17ed902 100644 --- a/src/sparse/sparse.jl +++ b/src/sparse/sparse.jl @@ -12,9 +12,9 @@ include("sparsematrixcsc.jl") # include("sparsevectordok.jl") # Old Julia syntax, not used, skipping for now mutable struct SparseMatrixCOO{T<:Real} - I :: Vector{Int} - J :: Vector{Int} - V :: Vector{T} + I::Vector{Int} + J::Vector{Int} + V::Vector{T} end const SparseVectorCOO = SparseMatrixCOO @@ -120,10 +120,10 @@ Matrix(A) """ function add!(A::SparseMatrixCOO, dofs1::AbstractVector{Int}, dofs2::AbstractVector{Int}, data) n, m = length(dofs1), length(dofs2) - @assert length(data) == n*m + @assert length(data) == n * m k = 1 - for j=1:m - for i=1:n + for j = 1:m + for i = 1:n add!(A, dofs1[i], dofs2[j], data[k]) k += 1 end @@ -145,7 +145,7 @@ function add!(A::SparseMatrixCOO, dofs::Vector{Int}, data::Array{Float64}, dim:: error("Simulation stopped.") end append!(A.I, dofs) - append!(A.J, dim*ones(Int, length(dofs))) + append!(A.J, dim * ones(Int, length(dofs))) append!(A.V, vec(data)) end diff --git a/src/sparse/sparsevectordok.jl b/src/sparse/sparsevectordok.jl index 01be819..d95297f 100644 --- a/src/sparse/sparsevectordok.jl +++ b/src/sparse/sparsevectordok.jl @@ -2,11 +2,11 @@ # License is MIT: see https://github.com/JuliaFEM/FEMSparse.jl/blob/master/LICENSE mutable struct SparseVectorDOK{Tv,Ti<:Integer} <: AbstractSparseArray{Tv,Ti,1} - data :: Dict{Ti,Tv} + data::Dict{Ti,Tv} end function SparseVectorDOK() - return SparseVectorDOK(Dict{Int64, Float64}()) + return SparseVectorDOK(Dict{Int64,Float64}()) end function SparseVectorDOK{Tv,Ti<:Integer}(b::SparseVector{Tv,Ti}) @@ -28,7 +28,7 @@ end function add!{Tv,Ti<:Integer}(b::SparseVectorDOK{Tv,Ti}, dofs::Vector{Ti}, data::Vector{Tv}) @assert length(dofs) == length(data) z = Tv(0) - for i=1:length(dofs) + for i = 1:length(dofs) @inbounds b.data[dofs[i]] = Base.get(b.data, i, z) + data[i] end return nothing @@ -39,7 +39,7 @@ function get{Tv,Ti<:Integer}(b::SparseVectorDOK{Tv,Ti}, i::Ti) end function get!{Tv,Ti<:Integer}(b::SparseVectorDOK{Tv,Ti}, dofs::Vector{Ti}, data::Vector{Tv}) - for (i,j) in enumerate(dofs) + for (i, j) in enumerate(dofs) data[i] = get(b, i) end return nothing