diff --git a/test/test_elasticity_2d_linear_with_surface_load.jl b/test/test_elasticity_2d_linear_with_surface_load.jl index ce6bb09..f0019b9 100644 --- a/test/test_elasticity_2d_linear_with_surface_load.jl +++ b/test/test_elasticity_2d_linear_with_surface_load.jl @@ -8,15 +8,22 @@ using JuliaFEM.Test @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...) + block = Problem(Elasticity, "BLOCK", 2) + block.properties.formulation = :plane_stress + block.properties.finite_strain = false + + elements = aster_create_elements(mesh, :BLOCK, :QU4) + update!(elements, "youngs modulus", 288.0) + update!(elements, "poissons ratio", 1/3) + update!(elements, "displacement load 2", 576.0) + push!(block, elements...) + + traction = aster_create_elements(mesh, :TOP, :SE2) + update!(traction, "displacement traction force 2", 288.0) + push!(block, traction...) + # boundary conditions bc_sym = Problem(Dirichlet, "symmetry bc", 2, "displacement") bc_elements_left = aster_create_elements(mesh, :LEFT, :SE2) @@ -24,14 +31,15 @@ using JuliaFEM.Test 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) + push!(solver, block, 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] + f = 288.0 + g = 576.0 + E = 288.0 + nu = 1/3 + u3_expected = f/E*[-nu, 1] + g/(2*E)*[-nu, 1] + u3 = reshape(block.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 index 2067196..b3426ec 100644 --- a/test/test_elasticity_3d_linear_with_surface_load.jl +++ b/test/test_elasticity_3d_linear_with_surface_load.jl @@ -18,9 +18,10 @@ using JuliaFEM.Test element1 = Element(Hex8, [1, 2, 3, 4, 5, 6, 7, 8]) element2 = Element(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]) + update!([element1], "youngs modulus", 288.0) + update!([element1], "poissons ratio", 1/3) + update!([element2], "displacement traction force 3", 288.0) + update!([element1], "displacement load 3", 576.0) elasticity_problem = Problem(Elasticity, "solve continuum block", 3) elasticity_problem.properties.finite_strain = false @@ -31,9 +32,9 @@ using JuliaFEM.Test symxz = Element(Quad4, [1, 2, 6, 5]) symyz = Element(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 + symxz["displacement 2"] = 0.0 + symxy["displacement 3"] = 0.0 boundary_problem = Problem(Dirichlet, "symmetry boundary conditions", 3, "displacement") push!(boundary_problem, symxy, symxz, symyz) @@ -44,6 +45,6 @@ using JuliaFEM.Test 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] + u_expected = 2.0 * [-1/3, -1/3, 1.0] @test isapprox(disp, u_expected) end diff --git a/test/test_elasticity_3d_unit_block.jl b/test/test_elasticity_3d_unit_block.jl new file mode 100644 index 0000000..b79343f --- /dev/null +++ b/test/test_elasticity_3d_unit_block.jl @@ -0,0 +1,99 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md + +using JuliaFEM +using JuliaFEM.Preprocess +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) + + block = Problem(Elasticity, fn, 3) + block.properties.finite_strain = false + + elements = aster_create_elements(mesh, :BLOCK, vol) + update!(elements, "youngs modulus", 288.0) + update!(elements, "poissons ratio", 1/3) + update!(elements, "displacement load 3", 576.0) + push!(block, elements...) + + traction = aster_create_elements(mesh, :LOAD, sur) + update!(traction, "displacement traction force 3", 288.0) + push!(block, traction...) + + bc = Problem(Dirichlet, "symmetry boundary condition", 3, "displacement") +# bc.properties.formulation = :incremental + symyz = aster_create_elements(mesh, :SYMYZ, sur) + symxz = aster_create_elements(mesh, :SYMXZ, sur) + symxy = aster_create_elements(mesh, :SYMXY, sur) + update!(symyz, "displacement 1", 0.0) + update!(symxz, "displacement 2", 0.0) + update!(symxy, "displacement 3", 0.0) + push!(bc, symyz..., symxz..., symxy...) + return block, bc, elements, traction, symyz, symxz, symxy +end + +function calc_size(elements, dim) + A = 0.0 + for element in elements + Ael = 0.0 + size(element, 1) == dim || continue + for (w, xi) in get_integration_points(element) + detJ = element(xi, 0.0, Val{:detJ}) + Ael += w*detJ + end + for (i, X) in enumerate(element["geometry"](0.0)) + info("$i : $X") + end + info("Area / volume: $Ael") + A += Ael + end + 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 +=# + +@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) + V = calc_size(block.elements, 3) + info("volume of block: $V") + A = calc_size(bc.elements, 2) + info("area of boundary condition: $A") + At = calc_size(traction, 2) + info("area of load surface: $At") + @test isapprox(V, 1.0) + @test isapprox(At, 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(block.assembly.u) + nu = round(Int, length(block.assembly.u)/3) + u = reshape(block.assembly.u, 3, nu) + f = reshape(full(block.assembly.f), 3, nu) + dump(round(u', 5)) + dump(round(f', 5)) + info("max |u| = $max_u") + @test isapprox(max_u, 2.0) +end + diff --git a/test/test_elasticity_tet10_stiffness_matrix.jl b/test/test_elasticity_tet10_stiffness_matrix.jl new file mode 100644 index 0000000..4cd07dd --- /dev/null +++ b/test/test_elasticity_tet10_stiffness_matrix.jl @@ -0,0 +1,33 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md + +using JuliaFEM +using JuliaFEM.Preprocess +using JuliaFEM.Test + +@testset "test tet10 stiffness matrix" begin + el = Element(Tet10) + el["youngs modulus"] = 480.0 + el["poissons ratio"] = 1/3 + x1 = [2.0, 3.0, 4.0] + x2 = [6.0, 3.0, 2.0] + x3 = [2.0, 5.0, 1.0] + x4 = [4.0, 3.0, 6.0] + x5 = 0.5*(x1+x2) + x6 = 0.5*(x2+x3) + x7 = 0.5*(x3+x1) + x8 = 0.5*(x1+x4) + x9 = 0.5*(x2+x4) + x10 = 0.5*(x3+x4) + el["geometry"] = Vector{Float64}[x1, x2, x3, x4, x5, x6, x7, x8, x9, x10] + pr = Problem(Elasticity, "tet10", 3) + Kt, f = assemble(pr, el, 0.0) + eigs = real(eigvals(Kt)) + eigs_expected = [8809.45, 4936.01, 2880.56, 2491.66, 2004.85, + 1632.49, 1264.32, 1212.42, 817.905, + 745.755, 651.034, 517.441, 255.1, 210.955, + 195.832, 104.008, 72.7562, 64.4376, 53.8515, + 23.8417, 16.6354, 9.54682, 6.93361, 2.22099, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] + @test isapprox(eigs, eigs_expected; atol=1.0e-2) +end diff --git a/test/test_elasticity_tet4_stiffness_matrix.jl b/test/test_elasticity_tet4_stiffness_matrix.jl new file mode 100644 index 0000000..3194605 --- /dev/null +++ b/test/test_elasticity_tet4_stiffness_matrix.jl @@ -0,0 +1,33 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md + +using JuliaFEM +using JuliaFEM.Preprocess +using JuliaFEM.Test + +@testset "test tet4 stiffness matrix" begin + el = Element(Tet4) + el["youngs modulus"] = 96.0 + el["poissons ratio"] = 1/3 + x1 = [2.0, 3.0, 4.0] + x2 = [6.0, 3.0, 2.0] + x3 = [2.0, 5.0, 1.0] + x4 = [4.0, 3.0, 6.0] + el["geometry"] = Vector{Float64}[x1, x2, x3, x4] + pr = Problem(Elasticity, "tet4", 3) + Kt, f = assemble(pr, el, 0.0) + Kt_expected = [ + 149.0 108.0 24.0 -1.0 6.0 12.0 -54.0 -48.0 0.0 -94.0 -66.0 -36.0 + 108.0 344.0 54.0 -24.0 104.0 42.0 -24.0 -216.0 -12.0 -60.0 -232.0 -84.0 + 24.0 54.0 113.0 0.0 30.0 35.0 0.0 -24.0 -54.0 -24.0 -60.0 -94.0 + -1.0 -24.0 0.0 29.0 -18.0 -12.0 -18.0 24.0 0.0 -10.0 18.0 12.0 + 6.0 104.0 30.0 -18.0 44.0 18.0 12.0 -72.0 -12.0 0.0 -76.0 -36.0 + 12.0 42.0 35.0 -12.0 18.0 29.0 0.0 -24.0 -18.0 0.0 -36.0 -46.0 + -54.0 -24.0 0.0 -18.0 12.0 0.0 36.0 0.0 0.0 36.0 12.0 0.0 + -48.0 -216.0 -24.0 24.0 -72.0 -24.0 0.0 144.0 0.0 24.0 144.0 48.0 + 0.0 -12.0 -54.0 0.0 -12.0 -18.0 0.0 0.0 36.0 0.0 24.0 36.0 + -94.0 -60.0 -24.0 -10.0 0.0 0.0 36.0 24.0 0.0 68.0 36.0 24.0 + -66.0 -232.0 -60.0 18.0 -76.0 -36.0 12.0 144.0 24.0 36.0 164.0 72.0 + -36.0 -84.0 -94.0 12.0 -36.0 -46.0 0.0 48.0 36.0 24.0 72.0 104.0] + @test isapprox(Kt, Kt_expected) +end diff --git a/test/test_elasticity_tet4_volume_load.jl b/test/test_elasticity_tet4_volume_load.jl new file mode 100644 index 0000000..cd6bc9a --- /dev/null +++ b/test/test_elasticity_tet4_volume_load.jl @@ -0,0 +1,34 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md + +using JuliaFEM +using JuliaFEM.Preprocess +using JuliaFEM.Test + +@testset "test tet4 element under volume load" begin + X1 = [2.0, 3.0, 4.0] + X2 = [6.0, 3.0, 2.0] + X3 = [2.0, 5.0, 1.0] + X4 = [4.0, 3.0, 6.0] + e1 = Element(Tet4, [1, 2, 3, 4]) + e1["geometry"] = Node[X1, X2, X3, X4] + e1["youngs modulus"] = 96.0 + e1["poissons ratio"] = 1/3 + e1["displacement load 3"] = 784.0/110.0 + e2 = Element(Tri3, [1, 2, 3]) + e2["geometry"] = Node[X2, X1, X3] + e2["displacement 1"] = 0.0 + e2["displacement 2"] = 0.0 + e2["displacement 3"] = 0.0 + p1 = Problem(Elasticity, "tetra", 3) + p2 = Problem(Dirichlet, "bc", 3, "displacement") + p1.properties.finite_strain = false + push!(p1, e1) + push!(p2, e2) + s = Solver() + push!(s, p1, p2) + call(s) + u_4 = p1.assembly.u[10:end] + u_expected = [-3.0/220.0, -9.0/220.0, 1.0/10.0] + @test isapprox(u_4, u_expected) +end