From 3df97fdca05274dff3441bf859fec3907cb56c03 Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Thu, 11 Feb 2016 12:07:25 +0200 Subject: [PATCH] fixed tests --- src/elasticity.jl | 4 +- test/test_elasticity.jl | 116 --------------------- test/test_elasticity_surface_load.jl | 147 ++++++++++++++++++++++++--- 3 files changed, 135 insertions(+), 132 deletions(-) diff --git a/src/elasticity.jl b/src/elasticity.jl index 3d4af76..915dd88 100644 --- a/src/elasticity.jl +++ b/src/elasticity.jl @@ -23,7 +23,7 @@ function get_unknown_field_name(::Type{Elasticity}) end function get_formulation_type(problem::Problem{Elasticity}) - info("INCREMENTAL FORMULATION") + # we are solving residual and add increment to previous solution vector return :incremental end @@ -81,7 +81,7 @@ function assemble{El<:Union{Tri3,Tri6,Quad4}}(problem::Problem{Elasticity}, elem nu 1-nu 0 0 0 (1-2*nu)/2] else - error("unknown plane formulation: $(props.formulation)") + error("unknown 2d formulation: $(props.formulation)") end S = D*[GL[1,1]; GL[2,2]; 2*GL[1,2]] # PK2 stress tensor in voigt notation diff --git a/test/test_elasticity.jl b/test/test_elasticity.jl index 7ba22e0..8249df2 100644 --- a/test/test_elasticity.jl +++ b/test/test_elasticity.jl @@ -43,120 +43,4 @@ function test_elasticity_volume_load() end #test_elasticity_volume_load() - -function test_elasticity_surface_load() - N = Vector[[0.0, 0.0], [1.0, 0.0], [0.0, 1.0], [1.0, 1.0]] - - element1 = Quad4([1, 2, 4, 3]) - element1["geometry"] = Vector[N[1], N[2], N[4], N[3]] - element1["youngs modulus"] = 900.0 - element1["poissons ratio"] = 0.25 - element1["displacement"] = (0.0 => Vector{Float64}[[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]]) - - element2 = Seg2([3, 4]) - element2["geometry"] = Vector[N[3], N[4]] - element2["displacement traction force"] = Vector[[0.0, -100.0], [0.0, -100.0]] - element2["displacement"] = (0.0 => Vector{Float64}[[0.0, 0.0], [0.0, 0.0]]) - - #free_dofs = [3, 5, 6, 8] - free_dofs = [3, 6, 7, 8] - problem = PlaneStressElasticityProblem() - push!(problem, element1) - push!(problem, element2) - solve!(problem, free_dofs, 0.0; max_iterations=10) - disp = element1("displacement", [1.0, 1.0], 0.0) - info("displacement at tip: $disp") - # verified using Code Aster. - @test isapprox(disp, [3.17431158889468E-02, -1.38591518927826E-01]) -end - - -function test_continuum_elasticity_with_surface_load() - nodes = JuliaFEM.Preprocess.aster_parse_nodes(""" - 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 - """) - - function set_geometry!(element, nodes) - element["geometry"] = Vector{Float64}[nodes[i] for i in get_connectivity(element)] - end - element1 = Hex8([1, 2, 3, 4, 5, 6, 7, 8]) - set_geometry!(element1, nodes) -# element1["youngs modulus"] = 900.0 -# element1["poissons ratio"] = 0.25 - element1["youngs modulus"] = 900.0 - element1["poissons ratio"] = 0.25 - element1["displacement"] = (0.0 => Vector{Float64}[[0.0, 0.0, 0.0] for i=1:8]) - - element2 = Quad4([5, 6, 7, 8]) - set_geometry!(element2, nodes) - element2["displacement traction force"] = Vector{Float64}[[0.0, 0.0, -100.0] for i=1:4] - element2["displacement"] = (0.0 => Vector{Float64}[[0.0, 0.0, 0.0] for i=1:4]) - - problem = ElasticityProblem() - push!(problem, element1) - push!(problem, element2) - - free_dofs = zeros(Bool, 8, 3) - x = 1 - y = 2 - z = 3 - free_dofs[2, x] = true - free_dofs[3, [x, y]] = true - free_dofs[4, y] = true - free_dofs[5, z] = true - free_dofs[6, [x, z]] = true - free_dofs[7, [x, y, z]] = true - free_dofs[8, [y, z]] = true - free_dofs = find(vec(free_dofs')) - info("free dofs: $free_dofs") - - info("initial force vector") - ass = JuliaFEM.Core.assemble(problem, 0.0) - info(reshape(full(ass.force_vector), 3, 8)) - info("initial stiffness matrix") - dump(round(Int, full(ass.stiffness_matrix))[free_dofs, free_dofs]) - solve!(problem, free_dofs, 0.0; max_iterations=10) - -#= - dx = Quad4([1, 4, 8, 5]) - dx["displacement 1"] = 0.0 - dy = Quad4([1, 5, 6, 2]) - dy["displacement 2"] = 0.0 - dz = Quad4([1, 2, 3, 4]) - dz["displacement 3"] = 0.0 - bc = DirichletProblem("displacement", 3) - for el in [dx, dy, dz] - set_geometry!(el, nodes) - push!(bc, el) - end - solver = JuliaFEM.Core.DirectSolver() - push!(solver, problem) - push!(solver, bc) - solver.dump_matrices = true - solver.name = "3d_hex8" - solver(0.0) -=# - - disp = element1("displacement", [1.0, 1.0, 1.0], 0.0) - info("displacement at tip: $disp") - info("displacement on element: ") - for (i, d) in enumerate(element1("displacement", 0.0)) - @printf "%d % f % f % f\n" [i;d]... - end - # 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]) - #@test isapprox(disp, [2.80559539222183E-03, 2.80559539222183E-03, -1.13019918093242E-02]) -end -#test_continuum_elasticity_with_surface_load() - end diff --git a/test/test_elasticity_surface_load.jl b/test/test_elasticity_surface_load.jl index f8afcc7..9d0059a 100644 --- a/test/test_elasticity_surface_load.jl +++ b/test/test_elasticity_surface_load.jl @@ -3,8 +3,10 @@ using JuliaFEM.Test using JuliaFEM.Core: Node, update!, Quad4, Seg2, Problem, Elasticity, Solver, Dirichlet +using JuliaFEM.Preprocess: aster_parse_nodes + +@testset "test 2d linear elasticity with surface load" begin -function get_test_problem() nodes = Dict{Int64, Node}( 1 => [0.0, 0.0], 2 => [1.0, 0.0], @@ -16,7 +18,7 @@ function get_test_problem() f = -E/10.0 expected = f/E*[-nu, 1] - # field problem is plane stress linear elasticity + # field problem element1 = Quad4([1, 2, 3, 4]) element2 = Seg2([3, 4]) update!([element1, element2], "geometry", nodes) @@ -29,7 +31,7 @@ function get_test_problem() elasticity_problem.properties.formulation = :plane_stress push!(elasticity_problem, element1, element2) - # dirichlet boundary condition, symmetry + # boundary condition, displacement symmetry sym13 = Seg2([1, 2]) sym23 = Seg2([4, 1]) update!([sym13, sym23], "geometry", nodes) @@ -38,18 +40,9 @@ function get_test_problem() # type, name, dimension, unknown_field_name boundary_problem = Problem(Dirichlet, "symmetry boundaries", 2, "displacement") push!(boundary_problem, sym13, sym23) - return elasticity_problem, boundary_problem -end - -@testset "test 2d linear with surface load." begin - - E = 288.0 - nu = 1.0/3.0 - f = -E/10.0 - expected = f/E*[-nu, 1] - elasticity_problem, boundary_problem = get_test_problem() + solver = Solver("solve block problem") - #solver.is_linear_system = true + solver.is_linear_system = true # to get linear solution push!(solver, elasticity_problem) push!(solver, boundary_problem) call(solver) @@ -60,3 +53,129 @@ end 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 = JuliaFEM.Preprocess.aster_parse_nodes(""" + 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 + """) + + function set_geometry!(element, nodes) + element["geometry"] = Vector{Float64}[nodes[i] for i in get_connectivity(element)] + end + element1 = Hex8([1, 2, 3, 4, 5, 6, 7, 8]) + set_geometry!(element1, nodes) + element1["youngs modulus"] = 900.0 + element1["poissons ratio"] = 0.25 + element1["displacement"] = (0.0 => Vector{Float64}[[0.0, 0.0, 0.0] for i=1:8]) + + element2 = Quad4([5, 6, 7, 8]) + set_geometry!(element2, nodes) + element2["displacement traction force"] = Vector{Float64}[[0.0, 0.0, -100.0] for i=1:4] + element2["displacement"] = (0.0 => Vector{Float64}[[0.0, 0.0, 0.0] for i=1:4]) + + problem = ElasticityProblem() + push!(problem, element1) + push!(problem, element2) + + free_dofs = zeros(Bool, 8, 3) + x = 1 + y = 2 + z = 3 + free_dofs[2, x] = true + free_dofs[3, [x, y]] = true + free_dofs[4, y] = true + free_dofs[5, z] = true + free_dofs[6, [x, z]] = true + free_dofs[7, [x, y, z]] = true + free_dofs[8, [y, z]] = true + free_dofs = find(vec(free_dofs')) + info("free dofs: $free_dofs") + + info("initial force vector") + ass = JuliaFEM.Core.assemble(problem, 0.0) + info(reshape(full(ass.force_vector), 3, 8)) + info("initial stiffness matrix") + dump(round(Int, full(ass.stiffness_matrix))[free_dofs, free_dofs]) + solve!(problem, free_dofs, 0.0; max_iterations=10) + +#= + dx = Quad4([1, 4, 8, 5]) + dx["displacement 1"] = 0.0 + dy = Quad4([1, 5, 6, 2]) + dy["displacement 2"] = 0.0 + dz = Quad4([1, 2, 3, 4]) + dz["displacement 3"] = 0.0 + bc = DirichletProblem("displacement", 3) + for el in [dx, dy, dz] + set_geometry!(el, nodes) + push!(bc, el) + end + solver = JuliaFEM.Core.DirectSolver() + push!(solver, problem) + push!(solver, bc) + solver.dump_matrices = true + solver.name = "3d_hex8" + solver(0.0) +=# + + disp = element1("displacement", [1.0, 1.0, 1.0], 0.0) + info("displacement at tip: $disp") + info("displacement on element: ") + for (i, d) in enumerate(element1("displacement", 0.0)) + @printf "%d % f % f % f\n" [i;d]... + end + # 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]) + #@test isapprox(disp, [2.80559539222183E-03, 2.80559539222183E-03, -1.13019918093242E-02]) +end +