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

most important tests pass now

Browse Source
This commit is contained in:
Jukka Aho
2016-05-22 17:00:01 +03:00
parent d6dea85b9d
commit 24656eb243
9 changed files with 107 additions and 67 deletions
Download Patch File Download Diff File Expand all files Collapse all files
src/JuliaFEM.jl
+1 -2
View File
@@ -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
src/elasticity.jl
+24 -28
View File
@@ -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
src/elements.jl
+19 -13
View File
@@ -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
src/integrate.jl
+6
View File
@@ -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
test/test_elasticity_2d_nonlinear_with_surface_load.jl
+2 -3
View File
@@ -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)
test/test_elasticity_2d_plane_stress_stiffness_matrix.jl
+38
View File
@@ -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
test/test_elasticity_3d_linear_with_surface_load.jl
+6 -7
View File
@@ -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
test/test_elasticity_3d_nonlinear_with_surface_load.jl
+6 -6
View File
@@ -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
test/test_elasticity_continuum.jl
+5 -8
View File
@@ -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)