From fe46907ac76be99fb5fd3a06602a4736a0f09892 Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Wed, 8 Jun 2016 21:50:45 +0300 Subject: [PATCH] buckling test, fixed bug with tensor <-> vector order --- src/elasticity.jl | 50 +++++++++++-------- .../test_elasticity_tet10_stiffness_matrix.jl | 2 +- test/test_elasticity_tet4_stiffness_matrix.jl | 11 +++- test/test_elasticity_tetra.jl | 32 ++++++++++++ 4 files changed, 72 insertions(+), 23 deletions(-) diff --git a/src/elasticity.jl b/src/elasticity.jl index 440ae8d..3ddce32 100644 --- a/src/elasticity.jl +++ b/src/elasticity.jl @@ -29,12 +29,20 @@ function assemble!(assembly::Assembly, problem::Problem{Elasticity}, element::El gdofs = get_gdofs(problem, element) if problem.properties.formulation in [:plane_stress, :plane_strain] Kt, f = assemble(problem, element, time, Val{:plane}) + add!(assembly.K, gdofs, gdofs, Kt) + add!(assembly.f, gdofs, f) + return + elseif problem.properties.formulation in [:continuum_buckling] + Km, Kg = assemble(problem, element, time, Val{:continuum_buckling}) + add!(assembly.K, gdofs, gdofs, Km) + add!(assembly.Kg, gdofs, gdofs, Kg) + return else Kt, f = assemble(problem, element, time, Val{problem.properties.formulation}) + add!(assembly.K, gdofs, gdofs, Kt) + add!(assembly.f, gdofs, f) + return end - add!(assembly.K, gdofs, gdofs, Kt) - add!(assembly.f, gdofs, f) - return Kt, f end """ Elasticity equations for 2d cases. """ @@ -206,12 +214,12 @@ function assemble{El<:Union{Tet4, Tet10, Hex8}}(problem::Problem{Elasticity}, el BL[1, 3*(i-1)+1] = dN[1,i] BL[2, 3*(i-1)+2] = dN[2,i] BL[3, 3*(i-1)+3] = dN[3,i] - BL[4, 3*(i-1)+1] = dN[2,i] + dN[1,i] - BL[4, 3*(i-1)+2] = dN[2,i] + dN[1,i] - BL[5, 3*(i-1)+2] = dN[3,i] + dN[2,i] - BL[5, 3*(i-1)+3] = dN[3,i] + dN[2,i] - BL[6, 3*(i-1)+1] = dN[1,i] + dN[3,i] - BL[6, 3*(i-1)+3] = dN[1,i] + dN[3,i] + BL[4, 3*(i-1)+1] = dN[2,i] + BL[4, 3*(i-1)+2] = dN[1,i] + BL[5, 3*(i-1)+2] = dN[3,i] + BL[5, 3*(i-1)+3] = dN[2,i] + BL[6, 3*(i-1)+1] = dN[3,i] + BL[6, 3*(i-1)+3] = dN[1,i] end E = element("youngs modulus", ip, time) @@ -276,12 +284,12 @@ function assemble{El<:Union{Tet4, Tet10, Hex8}}(problem::Problem{Elasticity}, el BL[1, 3*(i-1)+1] = dN[1,i] BL[2, 3*(i-1)+2] = dN[2,i] BL[3, 3*(i-1)+3] = dN[3,i] - BL[4, 3*(i-1)+1] = dN[2,i] + dN[1,i] - BL[4, 3*(i-1)+2] = dN[2,i] + dN[1,i] - BL[5, 3*(i-1)+2] = dN[3,i] + dN[2,i] - BL[5, 3*(i-1)+3] = dN[3,i] + dN[2,i] - BL[6, 3*(i-1)+1] = dN[1,i] + dN[3,i] - BL[6, 3*(i-1)+3] = dN[1,i] + dN[3,i] + BL[4, 3*(i-1)+1] = dN[2,i] + BL[4, 3*(i-1)+2] = dN[1,i] + BL[5, 3*(i-1)+2] = dN[3,i] + BL[5, 3*(i-1)+3] = dN[2,i] + BL[6, 3*(i-1)+1] = dN[3,i] + BL[6, 3*(i-1)+3] = dN[1,i] end fill!(BNL, 0.0) @@ -314,9 +322,9 @@ function assemble{El<:Union{Tet4, Tet10, Hex8}}(problem::Problem{Elasticity}, el S3[1,1] = stress_vec[1] S3[2,2] = stress_vec[2] S3[3,3] = stress_vec[3] - S3[2,3] = S3[3,2] = stress_vec[4] - S3[1,3] = S3[3,1] = stress_vec[5] - S3[1,2] = S3[2,1] = stress_vec[6] + S3[1,2] = S3[2,1] = stress_vec[4] + S3[2,3] = S3[3,2] = stress_vec[5] + S3[1,3] = S3[3,1] = stress_vec[6] S3[4:6,4:6] = S3[7:9,7:9] = S3[1:3,1:3] Km += w*BL'*D*BL @@ -428,9 +436,9 @@ function assemble{El<:Union{Tet4, Tet10, Hex8}}(problem::Problem{Elasticity}, el S3[1,1] = stress_vec[1] S3[2,2] = stress_vec[2] S3[3,3] = stress_vec[3] - S3[2,3] = S3[3,2] = stress_vec[4] - S3[1,3] = S3[3,1] = stress_vec[5] - S3[1,2] = S3[2,1] = stress_vec[6] + S3[1,2] = S3[2,1] = stress_vec[4] + S3[2,3] = S3[3,2] = stress_vec[5] + S3[1,3] = S3[3,1] = stress_vec[6] S3[4:6,4:6] = S3[7:9,7:9] = S3[1:3,1:3] if props.finite_strain diff --git a/test/test_elasticity_tet10_stiffness_matrix.jl b/test/test_elasticity_tet10_stiffness_matrix.jl index 4cd07dd..c75c88a 100644 --- a/test/test_elasticity_tet10_stiffness_matrix.jl +++ b/test/test_elasticity_tet10_stiffness_matrix.jl @@ -21,7 +21,7 @@ using JuliaFEM.Test 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) + Kt, f = assemble(pr, el, 0.0, Val{:continuum_linear}) eigs = real(eigvals(Kt)) eigs_expected = [8809.45, 4936.01, 2880.56, 2491.66, 2004.85, 1632.49, 1264.32, 1212.42, 817.905, diff --git a/test/test_elasticity_tet4_stiffness_matrix.jl b/test/test_elasticity_tet4_stiffness_matrix.jl index 3194605..d69fa4c 100644 --- a/test/test_elasticity_tet4_stiffness_matrix.jl +++ b/test/test_elasticity_tet4_stiffness_matrix.jl @@ -15,7 +15,7 @@ using JuliaFEM.Test 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, f = assemble(pr, el, 0.0, Val{:continuum_linear}) 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 @@ -29,5 +29,14 @@ using JuliaFEM.Test -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] + if !isapprox(Kt, Kt_expected) + info("Test failed") + info("Kt_expected") + dump(Kt_expected) + info("Kt") + dump(Kt) + end + @test isapprox(Kt, Kt_expected) + Kt, f = assemble(pr, el, 0.0, Val{:continuum}) @test isapprox(Kt, Kt_expected) end diff --git a/test/test_elasticity_tetra.jl b/test/test_elasticity_tetra.jl index c4b92da..a0a60d9 100644 --- a/test/test_elasticity_tetra.jl +++ b/test/test_elasticity_tetra.jl @@ -64,3 +64,35 @@ end @test isapprox(u_4, u_expected) end +@testset "test tet4 + buckling" begin + X = Dict{Int, Vector{Float64}}( + 1 => [2.0, 3.0, 4.0], + 2 => [6.0, 3.0, 2.0], + 3 => [2.0, 5.0, 1.0], + 4 => [4.0, 3.0, 6.0]) + u = Dict{Int, Vector{Float64}}( + 1 => [0.0, 0.0, 0.0], + 2 => [0.0, 0.0, 0.0], + 3 => [0.0, 0.0, 0.0], + 4 => [-0.25, -0.25, -0.25]) + e1 = Element(Tet4, [1, 2, 3, 4]) + update!(e1, "geometry", X) + update!(e1, "displacement", u) + update!(e1, "youngs modulus", 96.0) + update!(e1, "poissons ratio", 1/3) + p1 = Problem(Elasticity, "tetra", 3) + p1.properties.formulation = :continuum_buckling + push!(p1, e1) + assemble!(p1, 0.0) + free_dofs = [10, 11, 12] + Km = sparse(p1.assembly.K)[free_dofs, free_dofs] + Kg = sparse(p1.assembly.Kg)[free_dofs, free_dofs] + dump(full(Km)) + dump(full(Kg)) + la = sort(eigs(Km, -Kg)[1]) + la_expected = [1.0, 4.0] + info("la = $la") + info("la_expected = $(la_expected)") + @test isapprox(la, la_expected) +end +