From 0ff12a9beea883738200cef82ce69342b1c87a3b Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Tue, 24 May 2016 19:05:25 +0300 Subject: [PATCH] med node ordering --- src/JuliaFEM.jl | 5 +- src/elements.jl | 4 ++ src/integrate.jl | 2 +- src/preprocess_aster_reader.jl | 31 +++++++--- test/test_elasticity_3d_unit_block.jl | 87 ++++++++++++++++++++------- 5 files changed, 96 insertions(+), 33 deletions(-) diff --git a/src/JuliaFEM.jl b/src/JuliaFEM.jl index 60c447e..37d29fe 100644 --- a/src/JuliaFEM.jl +++ b/src/JuliaFEM.jl @@ -22,7 +22,7 @@ export DCTI ### ELEMENTS ### include("elements.jl") # common element routines -export Node, AbstractElement, Element, update! +export Node, AbstractElement, Element, update!, get_connectivity include("lagrange_macro.jl") # Continuous Galerkin (Lagrange) elements generated using macro export Seg2, Tri3, Tri6, Quad4, Hex8, Tet4, Tet10 @@ -77,6 +77,9 @@ end module Postprocess include("xdmf.jl") +export xdmf_new_temporal_collection, xdmf_new_grid, + xdmf_new_mesh!, xdmf_new_nodal_field!, + xdmf_save_model, xdmf_new_model end """ JuliaFEM testing routines. """ diff --git a/src/elements.jl b/src/elements.jl index d9d6b3a..6c13080 100644 --- a/src/elements.jl +++ b/src/elements.jl @@ -23,6 +23,10 @@ function setindex!(element::Element, data, field_name::ASCIIString) element.fields[field_name] = Field(data) end +function call(element::Element, field_name::ASCIIString, time=0.0) + return element[field_name](time) +end + function call(element::Element, xi::Vector, time=0.0) get_basis(element, xi, time) end diff --git a/src/integrate.jl b/src/integrate.jl index 7482657..345498b 100644 --- a/src/integrate.jl +++ b/src/integrate.jl @@ -207,7 +207,7 @@ function get_integration_points(element::LinearElement; order=2) get_integration_points(element, Val{order}) end -function get_integration_points(element::QuadraticElement; order=3) +function get_integration_points(element::QuadraticElement; order=2) get_integration_points(element, Val{order}) end diff --git a/src/preprocess_aster_reader.jl b/src/preprocess_aster_reader.jl index 047467b..98afed5 100644 --- a/src/preprocess_aster_reader.jl +++ b/src/preprocess_aster_reader.jl @@ -6,11 +6,15 @@ using JuliaFEM function aster_create_elements(mesh, element_set, element_type=nothing; reverse_connectivity=false) elements = Element[] - mapping = Dict(:QU4 => Quad4, :TR3 => Tri3, :SE2 => Seg2, :HE8 => Hex8, :TE4 => Tet4) + mapping = Dict( + :SE2 => Seg2, + :TR3 => Tri3, + :TR6 => Tri6, + :QU4 => Quad4, + :HE8 => Hex8, + :TE4 => Tet4, + :T10 => Tet10) for (elid, (eltype, elset, elcon)) in mesh["connectivity"] - if !haskey(mapping, eltype) - error("aster_create_elements: unknown element mapping $eltype") - end elset == element_set || continue if !isa(element_type, Void) if isa(element_type, Tuple) @@ -24,6 +28,9 @@ function aster_create_elements(mesh, element_set, element_type=nothing; reverse_ if reverse_connectivity elcon = reverse(elcon) end + if !haskey(mapping, eltype) + error("aster_create_elements: unknown element mapping $eltype") + end element = Element(mapping[eltype], elcon) push!(elements, element) end @@ -290,10 +297,16 @@ end # hex8 nodes rotating cw first in yz plane then x+1 global const med_elmap = Dict{Symbol, Vector{Int}}( - :HE8 => [4, 8, 7, 3, 1, 5, 6, 2], - :TE4 => [2, 3, 4, 1], + :SE2 => [1, 2], + :SE3 => [1, 2, 3], + :TR3 => [1, 2, 3], + :TR6 => [1, 2, 3, 4, 5, 6], :QU4 => [1, 2, 3, 4], -# :SE2 => [2, 1] + :HE8 => [4, 8, 7, 3, 1, 5, 6, 2], # ..? + :TE4 => [3, 2, 1, 4], + :T10 => [3, 2, 1, 4, 6, 5, 7, 10, 9, 8] +# :T10 => [3, 4, 1, 2, 10, 8, 7, 6, 9, 5] +# :T10 => [5, 9, 6, 7, 8, 10, 2, 1, 4, 3] ) function get_connectivity(med::MEDFile, elsets, mesh_name) @@ -318,8 +331,8 @@ function get_connectivity(med::MEDFile, elsets, mesh_name) if haskey(med_elmap, eltype) elco = elco[med_elmap[eltype]] else - #warn("no element mapping info found for element type $eltype") - #warn("consider this as a warning: element may have french nodal ordering") + warn("no element mapping info found for element type $eltype") + warn("consider this as a warning: element may have french nodal ordering") end d[element_ids[i]] = (eltype, elset, elco) end diff --git a/test/test_elasticity_3d_unit_block.jl b/test/test_elasticity_3d_unit_block.jl index b79343f..4ad82d2 100644 --- a/test/test_elasticity_3d_unit_block.jl +++ b/test/test_elasticity_3d_unit_block.jl @@ -3,11 +3,13 @@ using JuliaFEM using JuliaFEM.Preprocess +using JuliaFEM.Postprocess using JuliaFEM.Test function get_model(fn, vol, sur) meshfile = Pkg.dir("JuliaFEM")*"/geometry/3d_blocks/BLOCK.med" mesh = parse_aster_med_file(meshfile, fn) + info(mesh) block = Problem(Elasticity, fn, 3) block.properties.finite_strain = false @@ -52,28 +54,54 @@ function calc_size(elements, dim) return A end -#= -@testset "test 3d block hex8" begin - block, bc = get_model("BLOCK_HEX8", :HE8, :QU4) - V = calc_size(block.elements, 3) - info("volume of block: $V") - A = calc_size(bc.elements, 2) - info("area of boundary condition: $A") - @test isapprox(V, 1.0) - @test isapprox(A, 3.0) - solver = Solver("solver block problem") - solver.is_linear_system = true - push!(solver, block, bc) - call(solver) - max_u = maximum(abs(block.assembly.u)) - info("max |u| = $max_u") - @test isapprox(max_u, 2.0) -end -=# +function xdmf_dump(all_elements, eltype, elsym, filename="/tmp/xdmf_result.xmf") + info("$(length(all_elements)) elements.") + xdoc, xmodel = xdmf_new_model() + coll = xdmf_new_temporal_collection(xmodel) + grid = xdmf_new_grid(coll; time=0.0) -@testset "test 3d block TET4" begin - block, bc, elements, traction, symyz, symxz, symxy = get_model("BLOCK_TET4", :TE4, :TR3) -# block, bc = get_model("BLOCK_TET10", :T10, :TR6) + Xg = Dict{Int64, Vector{Float64}}() + ug = Dict{Int64, Vector{Float64}}() + nids = Dict{Int64, Int64}() + for element in all_elements + conn = get_connectivity(element) + for (i, c) in enumerate(conn) + nids[c] = c + end + X = element("geometry", 0.0) + for (i, c) in enumerate(conn) + Xg[c] = X[i] + end + haskey(element, "displacement") || continue + u = element("displacement", 0.0) + for (i, c) in enumerate(conn) + ug[c] = u[i] + end + end + perm = sort(collect(keys(Xg))) + nodes = Vector{Float64}[Xg[i] for i in perm] + disp = Vector{Float64}[ug[i] for i in perm] + nids = Int[nids[i] for i in perm] + inids = Dict{Int64, Int64}() + for (i, nid) in enumerate(nids) + inids[nid] = i + end + elements = [] + for element in all_elements + isa(element, eltype) || continue + conn = get_connectivity(element) + nconn = [inids[i] for i in conn] + push!(elements, (elsym, nconn)) + end + + xdmf_new_mesh!(grid, nodes, elements) + xdmf_new_nodal_field!(grid, "displacement", disp) + xdmf_save_model(xdoc, filename) + info("model dumped to $filename") +end + +function calc_model(model, volume_element, surface_element) + block, bc, elements, traction, symyz, symxz, symxy = get_model(model, volume_element, surface_element) V = calc_size(block.elements, 3) info("volume of block: $V") A = calc_size(bc.elements, 2) @@ -94,6 +122,21 @@ end dump(round(u', 5)) dump(round(f', 5)) info("max |u| = $max_u") - @test isapprox(max_u, 2.0) + return block, u end +@testset "test 3d block TET4" begin + block, u = calc_model("BLOCK_TET4", :TE4, :TR3) + @test isapprox(maximum(u), 2.1329516539440205) +end + +@testset "test 3d block TET10" begin + block, u = calc_model("BLOCK_TET10", :T10, :TR6) + @test isapprox(maximum(u), 2.13656216413056) +end + +@testset "test 3d block HEX8" begin + block, u = calc_model("BLOCK_HEX8", :HE8, :QU4) +# xdmf_dump(block.elements, Element{Hex8}, :Hex8, "/tmp/BLOCK_HEX8.xmf") + @test isapprox(maximum(u), 2.0) +end