From e7f29b3c5af625a99dcf982ee134d7de03c93cb2 Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Sat, 13 Feb 2016 01:15:02 +0200 Subject: [PATCH] playing with an idea, should we have separate file for every test to make development easier...? --- ..._elasticity_2d_linear_with_surface_load.jl | 38 +++++ ...asticity_2d_nonlinear_with_surface_load.jl | 35 +++++ ..._elasticity_3d_linear_with_surface_load.jl | 50 +++++++ ...asticity_3d_nonlinear_with_surface_load.jl | 51 +++++++ test/test_elasticity_surface_load.jl | 139 ------------------ 5 files changed, 174 insertions(+), 139 deletions(-) create mode 100644 test/test_elasticity_2d_linear_with_surface_load.jl create mode 100644 test/test_elasticity_2d_nonlinear_with_surface_load.jl create mode 100644 test/test_elasticity_3d_linear_with_surface_load.jl create mode 100644 test/test_elasticity_3d_nonlinear_with_surface_load.jl diff --git a/test/test_elasticity_2d_linear_with_surface_load.jl b/test/test_elasticity_2d_linear_with_surface_load.jl new file mode 100644 index 0000000..9779baf --- /dev/null +++ b/test/test_elasticity_2d_linear_with_surface_load.jl @@ -0,0 +1,38 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md + +using JuliaFEM.Test + +using JuliaFEM.Preprocess: aster_create_elements, parse_aster_med_file +using JuliaFEM.Core: Problem, Elasticity, Dirichlet, Solver, update! + +@testset "test 2d linear elasticity with surface load" begin + meshfile = "/geometry/2d_block/BLOCK_1elem.med" + mesh = parse_aster_med_file(Pkg.dir("JuliaFEM")*meshfile) + # field problem + body = Problem(Elasticity, "BLOCK", 2) + body.properties.formulation = :plane_stress + body_elements = aster_create_elements(mesh, :BLOCK, :QU4) + update!(body_elements, "youngs modulus", 900.0) + update!(body_elements, "poissons ratio", 0.25) + trac_elements = aster_create_elements(mesh, :TOP, :SE2) + update!(trac_elements, "displacement traction force 2", -100.0) + push!(body, body_elements..., trac_elements...) + # boundary conditions + bc_sym = Problem(Dirichlet, "symmetry bc", 2, "displacement") + bc_elements_left = aster_create_elements(mesh, :LEFT, :SE2) + bc_elements_bottom = aster_create_elements(mesh, :BOTTOM, :SE2) + update!(bc_elements_left, "displacement 1", 0.0) + update!(bc_elements_bottom, "displacement 2", 0.0) + push!(bc_sym, bc_elements_left..., bc_elements_bottom...) + solver = Solver("solve block problem") + body.properties.finite_strain = false + push!(solver, body, bc_sym) + call(solver) + f = -100.0 + E = 900.0 + nu = 0.25 + u3_expected = f/E*[-nu, 1] + u3 = reshape(body.assembly.u, 2, 4)[:,3] + @test isapprox(u3, u3_expected) +end diff --git a/test/test_elasticity_2d_nonlinear_with_surface_load.jl b/test/test_elasticity_2d_nonlinear_with_surface_load.jl new file mode 100644 index 0000000..be77555 --- /dev/null +++ b/test/test_elasticity_2d_nonlinear_with_surface_load.jl @@ -0,0 +1,35 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md + +using JuliaFEM.Test + +using JuliaFEM.Preprocess: aster_create_elements, parse_aster_med_file +using JuliaFEM.Core: Problem, Elasticity, Dirichlet, Solver, update! + +@testset "test 2d nonlinear elasticity with surface load" begin + meshfile = "/geometry/2d_block/BLOCK_1elem.med" + mesh = parse_aster_med_file(Pkg.dir("JuliaFEM")*meshfile) + # field problem + body = Problem(Elasticity, "BLOCK", 2) + body.properties.formulation = :plane_stress + body_elements = aster_create_elements(mesh, :BLOCK, :QU4) + update!(body_elements, "youngs modulus", 900.0) + update!(body_elements, "poissons ratio", 0.25) + trac_elements = aster_create_elements(mesh, :TOP, :SE2) + update!(trac_elements, "displacement traction force 2", -100.0) + push!(body, body_elements..., trac_elements...) + # boundary conditions + bc_sym = Problem(Dirichlet, "symmetry bc", 2, "displacement") + bc_elements_left = aster_create_elements(mesh, :LEFT, :SE2) + bc_elements_bottom = aster_create_elements(mesh, :BOTTOM, :SE2) + update!(bc_elements_left, "displacement 1", 0.0) + update!(bc_elements_bottom, "displacement 2", 0.0) + push!(bc_sym, bc_elements_left..., bc_elements_bottom...) + solver = Solver("solve block problem") + push!(solver, body, bc_sym) + call(solver) + # result is verified using code aster + u3_expected = [3.17431158889468E-02, -1.38591518927826E-01] + u3 = reshape(body.assembly.u, 2, 4)[:,3] + @test isapprox(u3, u3_expected) +end diff --git a/test/test_elasticity_3d_linear_with_surface_load.jl b/test/test_elasticity_3d_linear_with_surface_load.jl new file mode 100644 index 0000000..da4a21c --- /dev/null +++ b/test/test_elasticity_3d_linear_with_surface_load.jl @@ -0,0 +1,50 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md + +using JuliaFEM.Test + +using JuliaFEM.Core: update!, Problem, Elasticity, Dirichlet, Solver, Node, Hex8, Quad4 + +@testset "test continuum 3d linear elasticity with surface load" begin + nodes = Dict{Int64, Node}( + 1 => [0.0, 0.0, 0.0], + 2 => [1.0, 0.0, 0.0], + 3 => [1.0, 1.0, 0.0], + 4 => [0.0, 1.0, 0.0], + 5 => [0.0, 0.0, 1.0], + 6 => [1.0, 0.0, 1.0], + 7 => [1.0, 1.0, 1.0], + 8 => [0.0, 1.0, 1.0]) + + element1 = Hex8([1, 2, 3, 4, 5, 6, 7, 8]) + element2 = Quad4([5, 6, 7, 8]) + update!([element1, element2], "geometry", nodes) + update!([element1], "youngs modulus", 900.0) + update!([element1], "poissons ratio", 0.25) + update!([element2], "displacement traction force", Vector{Float64}[[0.0, 0.0, -100.0] for i=1:4]) + + elasticity_problem = Problem(Elasticity, "solve continuum block", 3) + elasticity_problem.properties.finite_strain = false + push!(elasticity_problem, element1) + push!(elasticity_problem, element2) + + symxy = Quad4([1, 2, 3, 4]) + symxz = Quad4([1, 2, 6, 5]) + symyz = Quad4([1, 4, 8, 5]) + update!([symxy, symxz, symyz], "geometry", nodes) + symxy["displacement 3"] = 0.0 + symxz["displacement 2"] = 0.0 + symyz["displacement 1"] = 0.0 + boundary_problem = Problem(Dirichlet, "symmetry boundary conditions", 3, "displacement") + push!(boundary_problem, symxy, symxz, symyz) + + solver = Solver("solve 3d block") + push!(solver, elasticity_problem) + push!(solver, boundary_problem) + call(solver) + + disp = element1("displacement", [1.0, 1.0, 1.0], 0.0) + info("displacement at tip: $disp") + u_expected = -100.0/900.0 * [-0.25, -0.25, 1.0] + @test isapprox(disp, u_expected) +end diff --git a/test/test_elasticity_3d_nonlinear_with_surface_load.jl b/test/test_elasticity_3d_nonlinear_with_surface_load.jl new file mode 100644 index 0000000..94cd5fe --- /dev/null +++ b/test/test_elasticity_3d_nonlinear_with_surface_load.jl @@ -0,0 +1,51 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md + +using JuliaFEM.Test +using JuliaFEM.Core: Node, update!, Quad4, Hex8, Problem, Elasticity, Solver, Dirichlet + +@testset "test continuum nonlinear elasticity with surface load" begin + + nodes = Dict{Int64, Node}( + 1 => [0.0, 0.0, 0.0], + 2 => [1.0, 0.0, 0.0], + 3 => [1.0, 1.0, 0.0], + 4 => [0.0, 1.0, 0.0], + 5 => [0.0, 0.0, 1.0], + 6 => [1.0, 0.0, 1.0], + 7 => [1.0, 1.0, 1.0], + 8 => [0.0, 1.0, 1.0]) + + element1 = Hex8([1, 2, 3, 4, 5, 6, 7, 8]) + element2 = Quad4([5, 6, 7, 8]) + update!([element1, element2], "geometry", nodes) + update!([element1], "youngs modulus", 900.0) + update!([element1], "poissons ratio", 0.25) + update!([element2], "displacement traction force", Vector{Float64}[[0.0, 0.0, -100.0] for i=1:4]) + + elasticity_problem = Problem(Elasticity, "solve continuum block", 3) + push!(elasticity_problem, element1) + push!(elasticity_problem, element2) + + symxy = Quad4([1, 2, 3, 4]) + symxz = Quad4([1, 2, 6, 5]) + symyz = Quad4([1, 4, 8, 5]) + update!([symxy, symxz, symyz], "geometry", nodes) + symxy["displacement 3"] = 0.0 + symxz["displacement 2"] = 0.0 + symyz["displacement 1"] = 0.0 + boundary_problem = Problem(Dirichlet, "symmetry boundary conditions", 3, "displacement") + push!(boundary_problem, symxy, symxz, symyz) + + solver = Solver("solve 3d block") + push!(solver, elasticity_problem) + push!(solver, boundary_problem) + call(solver) + + disp = element1("displacement", [1.0, 1.0, 1.0], 0.0) + info("displacement at tip: $disp") + # verified using Code Aster. + # 2015-12-12-continuum-elasticity/vim c3d_grot_gdep_traction_force.comm + @test isapprox(disp, [3.17431158889468E-02, 3.17431158889468E-02, -1.38591518927826E-01]) +end + diff --git a/test/test_elasticity_surface_load.jl b/test/test_elasticity_surface_load.jl index 959379c..e69de29 100644 --- a/test/test_elasticity_surface_load.jl +++ b/test/test_elasticity_surface_load.jl @@ -1,139 +0,0 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md - -using JuliaFEM.Test -using JuliaFEM.Core: Node, update!, Quad4, Seg2, Hex8, Problem, Elasticity, Solver, Dirichlet -using JuliaFEM.Preprocess: aster_parse_nodes - -@testset "test 2d linear elasticity with surface load" begin - - nodes = Dict{Int64, Node}( - 1 => [0.0, 0.0], - 2 => [1.0, 0.0], - 3 => [1.0, 1.0], - 4 => [0.0, 1.0]) - - E = 288.0 - nu = 1.0/3.0 - f = -E/10.0 - expected = f/E*[-nu, 1] - - # field problem - element1 = Quad4([1, 2, 3, 4]) - element2 = Seg2([3, 4]) - update!([element1, element2], "geometry", nodes) - update!(element1, "youngs modulus", E) - update!(element1, "poissons ratio", nu) -# update!(element2, "displacement traction force", [0.0, f]) - update!(element2, "displacement traction force", Vector{Float64}[[0.0, f], [0.0, f]]) - # type, name, dimension - elasticity_problem = Problem(Elasticity, "block", 2) - elasticity_problem.properties.formulation = :plane_stress - push!(elasticity_problem, element1, element2) - - # boundary condition, displacement symmetry - sym13 = Seg2([1, 2]) - sym23 = Seg2([4, 1]) - update!([sym13, sym23], "geometry", nodes) - update!(sym13, "displacement 2", 0.0) - update!(sym23, "displacement 1", 0.0) - # type, name, dimension, unknown_field_name - boundary_problem = Problem(Dirichlet, "symmetry boundaries", 2, "displacement") - push!(boundary_problem, sym13, sym23) - - solver = Solver("solve block problem") - solver.is_linear_system = true # to get linear solution - push!(solver, elasticity_problem) - push!(solver, boundary_problem) - call(solver) - element1 = elasticity_problem.elements[1] - u_disp = element1("displacement", [1.0, 1.0], 0.0) - info("Displacement = $u_disp") - @test isapprox(u_disp, expected) - -end - -@testset "test 2d nonlinear elasticity with surface load" begin - nodes = Dict{Int64, Node}( - 1 => [0.0, 0.0], - 2 => [1.0, 0.0], - 3 => [0.0, 1.0], - 4 => [1.0, 1.0]) - element1 = Quad4([1, 2, 4, 3]) - update!(element1, "geometry", nodes) - element1["youngs modulus"] = 900.0 - element1["poissons ratio"] = 0.25 - element2 = Seg2([3, 4]) - update!(element2, "geometry", nodes) - element2["displacement traction force"] = Vector[[0.0, -100.0], [0.0, -100.0]] - elasticity_problem = Problem(Elasticity, "bock", 2) - elasticity_problem.properties.formulation = :plane_stress - push!(elasticity_problem, element1, element2) - - # boundary condition, displacement symmetry - sym13 = Seg2([1, 2]) - sym23 = Seg2([3, 1]) - update!([sym13, sym23], "geometry", nodes) - update!(sym13, "displacement 2", 0.0) - update!(sym23, "displacement 1", 0.0) - # type, name, dimension, unknown_field_name - boundary_problem = Problem(Dirichlet, "symmetry boundaries", 2, "displacement") - push!(boundary_problem, sym13, sym23) - - solver = Solver("solve block problem") - push!(solver, elasticity_problem) - push!(solver, boundary_problem) - call(solver) - element1 = elasticity_problem.elements[1] - - u_disp = element1("displacement", [1.0, 1.0], 0.0) - # verified using Code Aster. - u_expected = [3.17431158889468E-02, -1.38591518927826E-01] - info("Displacement = $u_disp") - @test isapprox(u_disp, u_expected) -end - -@testset "test continuum linear elasticity with surface load" begin - - nodes = Dict{Int64, Node}( - 1 => [0.0, 0.0, 0.0], - 2 => [1.0, 0.0, 0.0], - 3 => [1.0, 1.0, 0.0], - 4 => [0.0, 1.0, 0.0], - 5 => [0.0, 0.0, 1.0], - 6 => [1.0, 0.0, 1.0], - 7 => [1.0, 1.0, 1.0], - 8 => [0.0, 1.0, 1.0]) - - element1 = Hex8([1, 2, 3, 4, 5, 6, 7, 8]) - element2 = Quad4([5, 6, 7, 8]) - update!([element1, element2], "geometry", nodes) - update!([element1], "youngs modulus", 900.0) - update!([element1], "poissons ratio", 0.25) - update!([element2], "displacement traction force", Vector{Float64}[[0.0, 0.0, -100.0] for i=1:4]) - - elasticity_problem = Problem(Elasticity, "solve continuum block", 3) - push!(elasticity_problem, element1) - push!(elasticity_problem, element2) - - symxy = Quad4([1, 2, 3, 4]) - symxz = Quad4([1, 2, 6, 5]) - symyz = Quad4([1, 4, 8, 5]) - update!([symxy, symxz, symyz], "geometry", nodes) - symxy["displacement 3"] = 0.0 - symxz["displacement 2"] = 0.0 - symyz["displacement 1"] = 0.0 - boundary_problem = Problem(Dirichlet, "symmetry boundary conditions", 3, "displacement") - push!(boundary_problem, symxy, symxz, symyz) - - solver = Solver("solve 3d block") - push!(solver, elasticity_problem) - push!(solver, boundary_problem) - call(solver) - - disp = element1("displacement", [1.0, 1.0, 1.0], 0.0) - info("displacement at tip: $disp") - # verified using Code Aster. - # 2015-12-12-continuum-elasticity/vim c3d_grot_gdep_traction_force.comm - @test isapprox(disp, [3.17431158889468E-02, 3.17431158889468E-02, -1.38591518927826E-01]) -end