diff --git a/.travis.yml b/.travis.yml index 0751371..2af6138 100644 --- a/.travis.yml +++ b/.travis.yml @@ -3,6 +3,10 @@ os: - linux julia: - 0.6 +addons: + apt: + packages: + - hdf5-tools notifications: email: false webhooks: @@ -17,4 +21,4 @@ before_script: script: - julia --color=yes -e 'using PkgTestSuite; test()' after_success: - - julia --color=yes -e 'using PkgTestSuite; deploy()' \ No newline at end of file + - julia --color=yes -e 'using PkgTestSuite; deploy()' diff --git a/deps/build.jl b/deps/build.jl index 97363c5..ee4c2f9 100644 --- a/deps/build.jl +++ b/deps/build.jl @@ -3,5 +3,7 @@ Pkg.clone("https://github.com/JuliaFEM/AbaqusReader.jl.git") Pkg.build("AbaqusReader") +Pkg.clone("https://github.com/JuliaFEM/AsterReader.jl.git") +Pkg.build("AsterReader") Pkg.clone("https://github.com/JuliaFEM/FEMQuad.jl.git") Pkg.build("FEMQuad") diff --git a/src/JuliaFEM.jl b/src/JuliaFEM.jl index 7649e38..11f1964 100644 --- a/src/JuliaFEM.jl +++ b/src/JuliaFEM.jl @@ -136,17 +136,13 @@ include("preprocess.jl") export create_elements, Mesh, add_node!, add_nodes!, add_element!, add_elements!, add_element_to_element_set!, add_node_to_node_set!, find_nearest_nodes, find_nearest_node, reorder_element_connectivity!, - create_node_set_from_element_set! + create_node_set_from_element_set!, filter_by_element_set include("preprocess_abaqus_reader.jl") export abaqus_read_mesh, create_surface_elements include("preprocess_aster_reader.jl") -export aster_create_elements, parse_aster_med_file, is_aster_mail_keyword, - parse_aster_header, aster_parse_nodes, aster_renumber_nodes!, - aster_renumber_elements!, aster_combine_meshes, aster_read_mesh, - filter_by_element_set, filter_by_element_id, MEDFile, aster_read_data, - aster_read_mesh_names, aster_read_node_sets, aster_read_nodes, RMEDFile +export aster_read_mesh end diff --git a/src/preprocess_aster_reader.jl b/src/preprocess_aster_reader.jl index ae9f91c..e2fb61e 100644 --- a/src/preprocess_aster_reader.jl +++ b/src/preprocess_aster_reader.jl @@ -1,206 +1,19 @@ # This file is a part of JuliaFEM. # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md -using HDF5 -using JuliaFEM - -function aster_parse_nodes(section; strip_characters=true) - nodes = Dict{Any, Vector{Float64}}() - has_started = false - for line in split(section, '\n') - m = matchall(r"[\w.-]+", line) - if (length(m) != 1) && (!has_started) - continue - end - 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 length(m) == 4 - nid = m[1] - if strip_characters - nid = matchall(r"\d", nid) - nid = parse(Int, nid[1]) - end - nodes[nid] = float(m[2:end]) - end - end - return nodes -end - - -""" Code Aster binary file (.med). """ -type 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. +using AsterReader """ -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 => ["NALL"]) - 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 +Nodes are different order in Code Aster compared to ABAQUS. This is yet +incomplete mapping between the permutations. Most of this is still a +great mystery. -""" Return element sets from med file. - -Notes ------ -One element set id can have multiple names. +# References +- http://onelab.info/pipermail/gmsh/2008/003850.html +- http://caelinux.org/wiki/index.php/Proj:UNVConvert """ -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 - for (k, v) in mesh["ELEME"] - elset_id = parse(Int, split(k, '_')[2]) - element_sets[elset_id] = collect(to_ascii(d) for d in v["GRO"]["NOM"]) - 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"] - node_ids = nodes["NUM"] - nset_ids = nodes["FAM"] - nnodes = length(node_ids) - 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"] - element_ids = elements["NUM"] - nelements = length(element_ids) - 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 - -""" Parse code aster .med file. - -Paramters ---------- -fn - file name to parse -mesh_name :: optional - mesh name, if several meshes in one file - -Returns -------- -Dict containing fields "nodes" and "connectivity". - -""" -function parse_aster_med_file(fn, mesh_name=nothing) - med = MEDFile(fn) - mesh_names = get_mesh_names(med::MEDFile) - 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 - - debug("Code Aster .med reader info:") - elsets = get_element_sets(med, mesh_name) - for (k,v) in elsets - debug("ELSET $k => $v") - end - nsets = get_node_sets(med, mesh_name) - for (k,v) in nsets - debug("NSET $k => $v") - end - nodes = get_nodes(med, nsets, mesh_name) - conn = get_connectivity(med, elsets, mesh_name) - result = Dict("nodes" => nodes, "connectivity" => conn) - result["nodes"] = nodes - result["connectivity"] = conn - return result -end - -# some glues about ordering, this is still a mystery.. -# http://onelab.info/pipermail/gmsh/2008/003850.html -# http://caelinux.org/wiki/index.php/Proj:UNVConvert - -#global const med_connectivity = Dict{Symbol, Vector{Int}}( -# :Tet4 => [3, 2, 1, 4], -# :Hex8 => [4, 8, 7, 3, 1, 5, 6, 2], # ..? -# :Tet10 => [3, 2, 1, 4, 6, 5, 7, 10, 9, 8]) - -global const med_connectivity = Dict{Symbol, Vector{Int}}( +const med_connectivity = Dict{Symbol, Vector{Int}}( :Tet4 => [4,3,1,2], :Tet10 => [4,3,1,2,10,7,8,9,6,5], :Pyr5 => [1,4,3,2,5], @@ -209,119 +22,46 @@ global const med_connectivity = Dict{Symbol, Vector{Int}}( :Hex20 => [4,8,7,3,1,5,6,2,20,15,19,11,12,16,14,10,17,13,18,9], :Hex27 => [4,8,7,3,1,5,6,2,20,15,19,11,12,16,14,10,17,13,18,9,24,25,26,23,21,22,27]) -# element names in CA -> element names in JuliaFEM -global const mapping = Dict( - +""" +Map element names used in in Code Aster to element names used in in JuliaFEM +""" +const med_element_names = Dict{Symbol, Symbol}( :PO1 => :Poi1, - :SE2 => :Seg2, :SE3 => :Seg3, :SE4 => :Seg4, - :TR3 => :Tri3, :TR6 => :Tri6, :TR7 => :Tri7, - :QU4 => :Quad4, :QU8 => :Quad8, :QU9 => :Quad9, - :TE4 => :Tet4, :T10 => :Tet10, - :PE6 => :Wedge6, :P15 => :Wedge15, :P18 => :Wedge18, - :HE8 => :Hex8, :H20 => :Hex20, :H27 => :Hex27, - :PY5 => :Pyr5, - :P13 => :Pyr13, + :P13 => :Pyr13) - ) +""" + aster_read_mesh(filename::String, mesh_name=nothing; reorder_element_connectivity=true) -""" Read code aster mesh and return Mesh. """ -function aster_read_mesh(fn, mesh_name=nothing; reorder_element_connectivity=true) - result = parse_aster_med_file(fn, mesh_name) - mesh = Mesh() - for (nid, (nsets, ncoords)) in result["nodes"] - add_node!(mesh, nid, ncoords) - for nset in nsets - add_node_to_node_set!(mesh, Symbol(nset), nid) - end - end - for (elid, (eltype, elsets, elcon)) in result["connectivity"] - haskey(mapping, eltype) || error("Code Aster .med reader: element type $eltype not found from mapping") - add_element!(mesh, elid, mapping[eltype], elcon) - for elset in elsets - add_element_to_element_set!(mesh, Symbol(elset), elid) - end +Read code aster mesh from file and return Mesh instance. If mesh file contains +several meshes, a name of mesh must be given. By default elements are reordered +so that they match to the conventions used in JuliaFEM. +""" +function aster_read_mesh(filename::String, mesh_name=nothing; reorder_element_connectivity=true) + m = AsterReader.aster_read_mesh(filename, mesh_name) + mesh = Mesh(m) + for (elid, eltype) in mesh.element_types + mesh.element_types[elid] = med_element_names[eltype] end if reorder_element_connectivity reorder_element_connectivity!(mesh, med_connectivity) end return mesh end - -""" Code Aster result file (.rmed). """ -type 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) - # INFO: quite safe assumption is that id is in node name, i.e. N1 => 1, N123 => 123 - node_id(node_name) = parse(matchall(r"\d+", node_name)[1]) - node_ids = map(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 contains(field_name, "ELGA") - 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/test/test_heat_3d_two_rings.jl b/test/test_heat_3d_two_rings.jl index 17ab52e..0d74776 100644 --- a/test/test_heat_3d_two_rings.jl +++ b/test/test_heat_3d_two_rings.jl @@ -5,6 +5,7 @@ using JuliaFEM using JuliaFEM.Preprocess using JuliaFEM.Postprocess using JuliaFEM.Testing +using AsterReader: RMEDFile, aster_read_nodes, aster_read_data #= Two rings, RING1 = inner, RING2 = outer, RINGS combined mesh. Set T=1.0 for diff --git a/test/test_mortar_3d_mesh_tie_two_rings.jl b/test/test_mortar_3d_mesh_tie_two_rings.jl index d6c1172..06e6067 100644 --- a/test/test_mortar_3d_mesh_tie_two_rings.jl +++ b/test/test_mortar_3d_mesh_tie_two_rings.jl @@ -6,6 +6,8 @@ using JuliaFEM.Preprocess using JuliaFEM.Postprocess using JuliaFEM.Testing +using AsterReader: RMEDFile, aster_read_nodes, aster_read_data + #= Two rings, RING1 is inner, RING2 is outer. Inner diameter is from 0.8 .. 0.9 and outer ring is 0.9 .. 1.0. Contact surface pair is RING1_OUTER <- RING2_INNER. diff --git a/test/test_preprocess.jl b/test/test_preprocess.jl index b2ad271..bb01a62 100644 --- a/test/test_preprocess.jl +++ b/test/test_preprocess.jl @@ -45,46 +45,6 @@ end @test first(nid) == 13 end -@testset "code aster / parse nodes" begin - section = """ - N9 2.0 3.0 4.0 - COOR_3D - N1 0.0 0.0 0.0 - N2 1.0 0.0 0.0 - N3 1.0 1.0 0.0 - N4 0.0 1.0 0.0 - N5 0.0 0.0 1.0 - N6 1.0 0.0 1.0 - N7 1.0 1.0 1.0 - N8 0.0 1.0 1.0 - FINSF - absdflasdf - N12 3.0 4.0 5.0 6.0 - N13 3.0 4.0 5.0 - """ - nodes = aster_parse_nodes(section) - @test nodes[1] == Float64[0.0, 0.0, 0.0] - @test nodes[8] == Float64[0.0, 1.0, 1.0] - @test length(nodes) == 8 -end - -@testset "test reading aster .med file" begin - meshfile = joinpath(datadir, "block_2d_1elem_quad4.med") - mesh = aster_read_mesh(meshfile) - @test length(mesh.element_sets) == 5 - @test length(mesh.node_sets) == 4 - @test length(mesh.elements) == 5 - @test length(mesh.nodes) == 4 - for elset in [:BLOCK, :TOP, :BOTTOM, :LEFT, :RIGHT] - @test haskey(mesh.element_sets, elset) - @test length(mesh.element_sets[elset]) == 1 - end - for nset in [:TOP_LEFT, :TOP_RIGHT, :BOTTOM_LEFT, :BOTTOM_RIGHT] - @test haskey(mesh.node_sets, nset) - @test length(mesh.node_sets[nset]) == 1 - end -end - @testset "test filter by element set" begin mesh = aster_read_mesh(joinpath(datadir, "block_2d_1elem_quad4.med")) mesh2 = filter_by_element_set(mesh, :BLOCK) @@ -127,30 +87,3 @@ end # @test isapprox(calculate_volume("PYRAMID_PYRAMID5_1", :Pyramid5, ?)) # @test isapprox(calculate_volume("PYRAMID_PYRAMID13_1", :Pyramid13, ?)) end - -@testset "read nodal field from code aster result file" begin - rmedfile = joinpath(datadir, "rings.rmed") - rmed = JuliaFEM.Preprocess.RMEDFile(rmedfile) - temp = JuliaFEM.Preprocess.aster_read_data(rmed, "TEMP") - @test isapprox(temp[15], 1.0) - @test isapprox(temp[95], 2.0) -end - -using JuliaFEM.Preprocess: MEDFile, get_element_sets - -@testset "test read element sets from med file, issue #111" begin - meshfile = joinpath(datadir, "hexmeshOverlappingGroups.med") - med = MEDFile(meshfile) - element_sets = get_element_sets(med, "Mesh_1") - @test element_sets[-10] == ["halfhex", "mosthex"] - @test element_sets[-11] == ["halfhex"] -end - -using JuliaFEM.Preprocess: aster_read_mesh - -@testset "test read overlapping ets, issue #111" begin - mesh_file = joinpath(datadir, "hexmeshOverlappingGroups.med") - mesh = aster_read_mesh(mesh_file, "Mesh_1") - @test length(mesh.element_sets[:mosthex]) == 273 - @test length(mesh.element_sets[:halfhex]) == 147 -end