From 24656eb24368ae24967e2d65a40e5de71f4c60e1 Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Sun, 22 May 2016 17:00:01 +0300 Subject: [PATCH] most important tests pass now --- src/JuliaFEM.jl | 3 +- src/elasticity.jl | 52 +++++++++---------- src/elements.jl | 32 +++++++----- src/integrate.jl | 6 +++ ...asticity_2d_nonlinear_with_surface_load.jl | 5 +- ...ticity_2d_plane_stress_stiffness_matrix.jl | 38 ++++++++++++++ ..._elasticity_3d_linear_with_surface_load.jl | 13 +++-- ...asticity_3d_nonlinear_with_surface_load.jl | 12 ++--- test/test_elasticity_continuum.jl | 13 ++--- 9 files changed, 107 insertions(+), 67 deletions(-) create mode 100644 test/test_elasticity_2d_plane_stress_stiffness_matrix.jl diff --git a/src/JuliaFEM.jl b/src/JuliaFEM.jl index 8f0a144..248daf3 100644 --- a/src/JuliaFEM.jl +++ b/src/JuliaFEM.jl @@ -13,7 +13,6 @@ autodiffcache = ForwardDiffCache() # export derivative, jacobian, hessian include("common.jl") -typealias Node Vector{Float64} include("fields.jl") export DCTI @@ -23,7 +22,7 @@ export DCTI ### ELEMENTS ### include("elements.jl") # common element routines -export Element, update! +export Node, Element, update! include("lagrange_macro.jl") # Continuous Galerkin (Lagrange) elements generated using macro export Seg2, Tri3, Quad4, Hex8, Tet4 diff --git a/src/elasticity.jl b/src/elasticity.jl index 4557022..c4459ea 100644 --- a/src/elasticity.jl +++ b/src/elasticity.jl @@ -64,8 +64,7 @@ function assemble{El<:Union{Tri3,Tri6,Quad4}}(problem::Problem{Elasticity}, elem for (w, xi) in get_integration_points(element) - J = element(xi, time, Val{:Jacobian}) - w = w*det(J) + detJ = element(xi, time, Val{:detJ}) N = element(xi, time) dN = element(xi, time, Val{:Grad}) @@ -126,16 +125,16 @@ function assemble{El<:Union{Tri3,Tri6,Quad4}}(problem::Problem{Elasticity}, elem S2[1,2] = S2[2,1] = S[3] S2[3:4,3:4] = S2[1:2,1:2] - Kt += w*BL'*D*BL # material stiffness + Kt += w*BL'*D*BL*detJ # material stiffness if props.finite_strain # add geometric stiffness - Kt += w*BNL'*S2*BNL # geometric stiffness + Kt += w*BNL'*S2*BNL*detJ # geometric stiffness end - f -= w*BL'*S # internal force + f -= w*BL'*S*detJ # internal force # volume load if haskey(element, "displacement load") b = element("displacement load", xi, time) - f += vec(w*N'*b) + f += w*vec(N'*b)*detJ end end @@ -153,8 +152,7 @@ function assemble{El<:Union{Seg2,Seg3}}(problem::Problem{Elasticity}, element::E for (w, xi) in get_integration_points(element) - J = element(xi, time, Val{:Jacobian}) - detJ = norm(J) + detJ = element(xi, time, Val{:detJ}) N = element(xi, time) if haskey(element, "displacement traction force") @@ -194,16 +192,15 @@ function assemble{El<:Union{Tet4, Tet10, Hex8}}(problem::Problem{Elasticity}, el Kt = zeros(dim*nnodes, dim*nnodes) f = zeros(dim*nnodes) - for ip in get_integration_points(element) - J = get_jacobian(element, ip, time) - w = ip.weight*det(J) - N = element(ip, time) - dN = element(ip, time, Val{:grad}) + for (w, xi) in get_integration_points(element) + detJ = element(xi, time, Val{:detJ}) + N = element(xi, time) + dN = element(xi, time, Val{:Grad}) # kinematics; calculate deformation gradient and strain gradu = zeros(dim, dim) if haskey(element, "displacement") - gradu += element("displacement", ip, time, Val{:grad}) + gradu += element("displacement", xi, time, Val{:Grad}) end strain = zeros(dim , dim) strain += 1/2*(gradu' + gradu) @@ -213,8 +210,8 @@ function assemble{El<:Union{Tet4, Tet10, Hex8}}(problem::Problem{Elasticity}, el strain += 1/2*gradu'*gradu end - E = element("youngs modulus", ip, time) - nu = element("poissons ratio", ip, time) + E = element("youngs modulus", xi, time) + nu = element("poissons ratio", xi, time) a = 1 - nu b = 1 - 2*nu @@ -273,16 +270,16 @@ function assemble{El<:Union{Tet4, Tet10, Hex8}}(problem::Problem{Elasticity}, el S3[1,2] = S3[2,1] = S[6] S3[4:6,4:6] = S3[7:9,7:9] = S3[1:3,1:3] - Kt += w*BL'*D*BL + Kt += w*BL'*D*BL*detJ if props.finite_strain - Kt += w*BNL'*S3*BNL + Kt += w*BNL'*S3*BNL*detJ end - f -= w*BL'*S + f -= w*BL'*S*detJ # volume load if haskey(element, "displacement load") T = element("displacement load", ip, time) - f += vec(w*T*N) + f += w*vec(T*N)*detJ end end @@ -298,18 +295,17 @@ function assemble{El<:Union{Tri3, Tri6, Quad4}}(problem::Problem{Elasticity}, el Kt = zeros(dim*nnodes, dim*nnodes) f = zeros(dim*nnodes) - for ip in get_integration_points(element) - JT = transpose(get_jacobian(element, ip, time)) - N = element(ip, time) - w = ip.weight*norm(cross(JT[:,1], JT[:,2])) + for (w, xi) in get_integration_points(element) + detJ = element(xi, time, Val{:detJ}) + N = element(xi, time) if haskey(element, "displacement traction force") - T = element("displacement traction force", ip, time) - f += vec(w*T*N) + T = element("displacement traction force", xi, time) + f += w*vec(T*N)*detJ end for i in 1:dim if haskey(element, "displacement traction force $i") - T = element("displacement traction force $i", ip, time) - f[i:dim:end] += vec(w*T*N) + T = element("displacement traction force $i", xi, time) + f[i:dim:end] += w*vec(T*N)*detJ end end end diff --git a/src/elements.jl b/src/elements.jl index 8d486ae..234eeb4 100644 --- a/src/elements.jl +++ b/src/elements.jl @@ -1,10 +1,10 @@ # This file is a part of JuliaFEM. # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md -import Base: getindex, setindex!, convert, size, length - abstract AbstractElement +typealias Node Vector{Float64} + type Element{E<:AbstractElement} connectivity :: Vector{Int} fields :: Dict{ASCIIString, Field} @@ -49,6 +49,20 @@ function call(element::Element, xi::Vector, time, ::Type{Val{:Jacobian}}) return J end +function call(element::Element, xi::Vector, time, ::Type{Val{:detJ}}) + J = element(xi, time, Val{:Jacobian}) + n, m = size(J) + if n == m # volume element + return det(J) + end + JT = transpose(J) + if size(JT, 2) == 1 # boundary of 2d problem, || ∂X/∂ξ || + return norm(JT) + else # manifold on 3d problem, || ∂X/∂ξ₁ × ∂X/∂ξ₂ || + return norm(cross(JT[:,1], JT[:,2])) + end +end + function get_jacobian{E}(element::Element{E}, xi::Vector, time=0.0) element(xi, time, Val{:Jacobian}) end @@ -128,18 +142,10 @@ function get_dualbasis(element::Element, time) De = zeros(nnodes, nnodes) Me = zeros(nnodes, nnodes) for (w, xi) in get_integration_points(element, Val{3}) - J = element(xi, time, Val{:Jacobian}) - JT = transpose(J) - if size(JT, 2) == 1 # plane problem - # || ∂X/∂ξ || - w *= norm(JT) - else - # || ∂X/∂ξ₁ × ∂X/∂ξ₂ || - w *= norm(cross(JT[:,1], JT[:,2])) - end + detJ = element(xi, time, Val{:detJ}) N = element(xi, time) - De += w*diagm(vec(N)) - Me += w*N'*N + De += w*diagm(vec(N))*detJ + Me += w*N'*N*detJ end return De, Me, De*inv(Me) end diff --git a/src/integrate.jl b/src/integrate.jl index fd8579a..08b187e 100644 --- a/src/integrate.jl +++ b/src/integrate.jl @@ -68,10 +68,12 @@ function get_integration_points(element::LineElement, ::Type{Val{3}}) [ (w[i], [xi[i]]) for i=1:3 ] end +#= function get_integration_points{E<:LineElement}(element::Element{E}, ::Type{Val{3}}) w, xi = get_integration_points(Val{3}) [ (w[i], [xi[i]]) for i=1:3 ] end +=# function get_integration_points(element::Quad4, ::Type{Val{2}}) w, xi = get_integration_points(Val{2}) @@ -106,6 +108,10 @@ function get_integration_points(element::Hex8) get_integration_points(element, Val{2}) end +function get_integration_points{E}(element::Element{E}, ::Type{Val{3}}) + get_integration_points(element.properties, Val{3}) +end + ### triangular and tetrahedral elements # http://math2.uncc.edu/~shaodeng/TEACHING/math5172/Lectures/Lect_15.PDF diff --git a/test/test_elasticity_2d_nonlinear_with_surface_load.jl b/test/test_elasticity_2d_nonlinear_with_surface_load.jl index be77555..2370740 100644 --- a/test/test_elasticity_2d_nonlinear_with_surface_load.jl +++ b/test/test_elasticity_2d_nonlinear_with_surface_load.jl @@ -1,11 +1,10 @@ # 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 -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) diff --git a/test/test_elasticity_2d_plane_stress_stiffness_matrix.jl b/test/test_elasticity_2d_plane_stress_stiffness_matrix.jl new file mode 100644 index 0000000..e42cd1d --- /dev/null +++ b/test/test_elasticity_2d_plane_stress_stiffness_matrix.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 + +# http://ahojukka5.github.io/posts/finite-element-solution-for-one-element-problem/ + +using JuliaFEM +using JuliaFEM.Test + +@testset "test 2d linear elasticity local matrices" begin + element = Element(Quad4) + element["geometry"] = Vector{Float64}[ + [0.0, 0.0], + [1.0, 0.0], + [1.0, 1.0], + [0.0, 1.0]] + element["youngs modulus"] = 288.0 + element["poissons ratio"] = 1/3 + element["displacement load"] = DCTI([4.0, 8.0]) + + problem = Problem(Elasticity, "[0x1] x [0x1] block", 2) + problem.properties.formulation = :plane_stress + K, f = assemble(problem, element) + + K_expected = [ + 144 54 -90 0 -72 -54 18 0 + 54 144 0 18 -54 -72 0 -90 + -90 0 144 -54 18 0 -72 54 + 0 18 -54 144 0 -90 54 -72 + -72 -54 18 0 144 54 -90 0 + -54 -72 0 -90 54 144 0 18 + 18 0 -72 54 -90 0 144 -54 + 0 -90 54 -72 0 18 -54 144] + + f_expected = [1, 2, 1, 2, 1, 2, 1, 2] + + @test isapprox(K, K_expected) + @test isapprox(f, f_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 da4a21c..2067196 100644 --- a/test/test_elasticity_3d_linear_with_surface_load.jl +++ b/test/test_elasticity_3d_linear_with_surface_load.jl @@ -1,10 +1,9 @@ # 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.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], @@ -16,8 +15,8 @@ using JuliaFEM.Core: update!, Problem, Elasticity, Dirichlet, Solver, Node, Hex8 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]) + 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) @@ -28,9 +27,9 @@ using JuliaFEM.Core: update!, Problem, Elasticity, Dirichlet, Solver, Node, Hex8 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]) + symxy = Element(Quad4, [1, 2, 3, 4]) + 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 diff --git a/test/test_elasticity_3d_nonlinear_with_surface_load.jl b/test/test_elasticity_3d_nonlinear_with_surface_load.jl index 94cd5fe..25349ad 100644 --- a/test/test_elasticity_3d_nonlinear_with_surface_load.jl +++ b/test/test_elasticity_3d_nonlinear_with_surface_load.jl @@ -1,8 +1,8 @@ # 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.Test -using JuliaFEM.Core: Node, update!, Quad4, Hex8, Problem, Elasticity, Solver, Dirichlet @testset "test continuum nonlinear elasticity with surface load" begin @@ -16,8 +16,8 @@ using JuliaFEM.Core: Node, update!, Quad4, Hex8, Problem, Elasticity, Solver, Di 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]) + 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) @@ -27,9 +27,9 @@ using JuliaFEM.Core: Node, update!, Quad4, Hex8, Problem, Elasticity, Solver, Di 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]) + symxy = Element(Quad4, [1, 2, 3, 4]) + 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 diff --git a/test/test_elasticity_continuum.jl b/test/test_elasticity_continuum.jl index 32bd07b..4e97f5f 100644 --- a/test/test_elasticity_continuum.jl +++ b/test/test_elasticity_continuum.jl @@ -4,9 +4,6 @@ using JuliaFEM using JuliaFEM.Test -using JuliaFEM.Core: Quad4, Hex8, Elasticity, Dirichlet, Problem, Node, Solver -using JuliaFEM.Core: update! - @testset "test simple continuum block with surface traction" begin nodes = Dict{Int64, Node}( @@ -18,8 +15,8 @@ using JuliaFEM.Core: update! 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]) + 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) @@ -56,11 +53,11 @@ using JuliaFEM.Core: update! dump(u) =# - dx = Quad4([1, 4, 8, 5]) + dx = Element(Quad4, [1, 4, 8, 5]) dx["displacement 1"] = 0.0 - dy = Quad4([1, 5, 6, 2]) + dy = Element(Quad4, [1, 5, 6, 2]) dy["displacement 2"] = 0.0 - dz = Quad4([1, 2, 3, 4]) + dz = Element(Quad4, [1, 2, 3, 4]) dz["displacement 3"] = 0.0 bc = Problem(Dirichlet, "symmetries", 3, "displacement") update!([dx, dy, dz], "geometry", nodes)