Glassvision/JuliaFEM.jl
1
0
Fork 0
mirror of https://github.com/JuliaFEM/JuliaFEM.jl.git synced 2026-09-17 09:12:09 +00:00
Code Issues Packages Projects Releases Wiki Activity

buckling test, fixed bug with tensor <-> vector order

Browse Source
This commit is contained in:
Jukka Aho
2016-06-08 21:50:45 +03:00
parent b65468c5d1
commit fe46907ac7
4 changed files with 72 additions and 23 deletions
Download Patch File Download Diff File Expand all files Collapse all files
src/elasticity.jl
+29 -21
View File
@@ -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
test/test_elasticity_tet10_stiffness_matrix.jl
+1 -1
View File
@@ -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,
test/test_elasticity_tet4_stiffness_matrix.jl
+10 -1
View File
@@ -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
test/test_elasticity_tetra.jl
+32
View File
@@ -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