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

fixed bug related to strain component ordering and internal force vector

Browse Source
This commit is contained in:
Jukka Aho
2016-05-27 17:45:12 +03:00
parent d606bde168
commit ccbd2d7224
5 changed files with 63 additions and 40 deletions
Download Patch File Download Diff File Expand all files Collapse all files
src/JuliaFEM.jl
+2 -1
View File
@@ -35,10 +35,11 @@ include("integrate.jl") # default integration points for elements
export get_integration_points
include("sparse.jl")
export add!
include("problems.jl") # common problem routines
export Problem, AbstractProblem, FieldProblem, BoundaryProblem,
get_unknown_field_dimension
get_unknown_field_dimension, get_gdofs
include("elasticity.jl") # elasticity equations
export Elasticity
src/dirichlet.jl
+2 -1
View File
@@ -7,11 +7,12 @@ using incremental formulation.
"""
type Dirichlet <: BoundaryProblem
formulation :: Symbol
variational :: Bool
dual_basis :: Bool
end
function Dirichlet()
Dirichlet(:total, true)
Dirichlet(:incremental, true, false)
end
function get_unknown_field_name(::Type{Dirichlet})
src/elasticity.jl
+42 -29
View File
@@ -86,7 +86,7 @@ function assemble{El<:Union{Tri3,Tri6,Quad4}}(problem::Problem{Elasticity}, elem
error("unknown plane formulation: $(props.formulation)")
end
# calculate stress
strain_vec = [strain[1,1]; strain[2,2]; 2*strain[1,2]]
strain_vec = [strain[1,1]; strain[2,2]; 2.0*strain[1,2]]
stress_vec = D*strain_vec
stress = [stress_vec[1] stress_vec[3]; stress_vec[3] stress_vec[2]]
cauchy_stress = F'*stress*F/det(F)
@@ -195,15 +195,16 @@ function assemble{El<:Union{Tet4, Tet10, Hex8}}(problem::Problem{Elasticity}, el
Kt = zeros(dim*nnodes, dim*nnodes)
f = zeros(dim*nnodes)
for (w, xi) in get_integration_points(element)
detJ = element(xi, time, Val{:detJ})
N = element(xi, time)
dN = element(xi, time, Val{:Grad})
for ip in get_integration_points(element)
detJ = element(ip, time, Val{:detJ})
w = ip.weight*detJ
N = element(ip, time)
dN = element(ip, time, Val{:Grad})
# kinematics; calculate deformation gradient and strain
gradu = zeros(dim, dim)
if haskey(element, "displacement")
gradu += element("displacement", xi, time, Val{:Grad})
gradu += element("displacement", ip, time, Val{:Grad})
end
strain = zeros(dim , dim)
strain += 1/2*(gradu' + gradu)
@@ -213,8 +214,8 @@ function assemble{El<:Union{Tet4, Tet10, Hex8}}(problem::Problem{Elasticity}, el
strain += 1/2*gradu'*gradu
end
E = element("youngs modulus", xi, time)
nu = element("poissons ratio", xi, time)
E = element("youngs modulus", ip, time)
nu = element("poissons ratio", ip, time)
D = E/((1.0+nu)*(1.0-2.0*nu)) * [
1.0-nu nu nu 0.0 0.0 0.0
@@ -224,8 +225,11 @@ function assemble{El<:Union{Tet4, Tet10, Hex8}}(problem::Problem{Elasticity}, el
0.0 0.0 0.0 0.0 0.5-nu 0.0
0.0 0.0 0.0 0.0 0.0 0.5-nu]
# # PK2 stress tensor in voigt notation
S = D*[strain[1,1]; strain[2,2]; strain[3,3]; 2*strain[2,3]; 2*strain[1,3]; 2*strain[1,2]]
# PK2 stress tensor in voigt notation
# strain_vec = [strain[1,1]; strain[2,2]; strain[3,3]; 2*strain[2,3]; 2*strain[1,3]; 2*strain[1,2]]
# order 11, 22, 33, 12, 23, 13 is in many text books ..?
strain_vec = [strain[1,1]; strain[2,2]; strain[3,3]; 2*strain[1,2]; 2*strain[2,3]; 2*strain[1,3]]
stress_vec = D*strain_vec
# add contributions: material and geometric stiffness + internal forces
fill!(BL, 0.0)
@@ -263,36 +267,44 @@ function assemble{El<:Union{Tet4, Tet10, Hex8}}(problem::Problem{Elasticity}, el
BNL[9, 3*(i-1)+3] = dN[3,i]
end
S3 = zeros(3*dim, 3*dim)
S3[1,1] = S[1]
S3[2,2] = S[2]
S3[3,3] = S[3]
S3[2,3] = S3[3,2] = S[4]
S3[1,3] = S3[3,1] = S[5]
S3[1,2] = S3[2,1] = S[6]
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[4:6,4:6] = S3[7:9,7:9] = S3[1:3,1:3]
Kt += w*BL'*D*BL*detJ
Kt += w*BL'*D*BL
if props.finite_strain
Kt += w*BNL'*S3*BNL*detJ
Kt += w*BNL'*S3*BNL
end
if get_formulation_type(problem) == :incremental
f -= w*BL'*S*detJ
f -= w*BL'*stress_vec
end
# volume load
if haskey(element, "displacement load")
T = element("displacement load", ip, time)
f += w*vec(T*N)*detJ
f += w*vec(T*N)
end
for i=1:dim
if haskey(element, "displacement load $i")
b = element("displacement load $i", xi, time)
f[i:dim:end] += w*vec(b*N)*detJ
b = element("displacement load $i", ip, time)
f[i:dim:end] += w*vec(b*N)
end
end
end
#=
if get_formulation_type(problem) == :incremental
if haskey(element, "displacement")
f -= Kt*vec(element["displacement"](time))
end
end
=#
return Kt, f
end
@@ -305,17 +317,18 @@ function assemble{El<:Union{Tri3, Tri6, Quad4}}(problem::Problem{Elasticity}, el
Kt = zeros(dim*nnodes, dim*nnodes)
f = zeros(dim*nnodes)
for (w, xi) in get_integration_points(element)
detJ = element(xi, time, Val{:detJ})
N = element(xi, time)
for ip in get_integration_points(element)
detJ = element(ip, time, Val{:detJ})
w = ip.weight*detJ
N = element(ip, time)
if haskey(element, "displacement traction force")
T = element("displacement traction force", xi, time)
f += w*vec(T*N)*detJ
T = element("displacement traction force", ip, time)
f += w*vec(T*N)
end
for i in 1:dim
if haskey(element, "displacement traction force $i")
T = element("displacement traction force $i", xi, time)
f[i:dim:end] += w*vec(T*N)*detJ
T = element("displacement traction force $i", ip, time)
f[i:dim:end] += w*vec(T*N)
end
end
end
src/elements.jl
+16 -8
View File
@@ -11,11 +11,10 @@ type Element{E<:AbstractElement}
properties :: E
end
function Element{E<:AbstractElement}(::Type{E}, connectivity=[], id=-1, fields=Dict(), properties...)
function Element{E<:AbstractElement}(::Type{E}, connectivity=[], integration_points=[], id=-1, fields=Dict(), properties...)
variant = E(properties...)
ips = get_integration_points(variant)
integration_points = [IP(i, w, xi) for (i, (w, xi)) in enumerate(ips)]
Element{E}(id, connectivity, integration_points, fields, variant)
element = Element{E}(id, connectivity, integration_points, fields, variant)
return element
end
function getindex(element::Element, field_name::ASCIIString)
@@ -26,11 +25,11 @@ function setindex!(element::Element, data, field_name::ASCIIString)
element.fields[field_name] = Field(data)
end
function call(element::Element, field_name::ASCIIString, time=0.0)
function call(element::Element, field_name::ASCIIString, time)
return element[field_name](time)
end
function call(element::Element, ip, time=0.0)
function call(element::Element, ip, time)
get_basis(element, ip, time)
end
@@ -64,8 +63,12 @@ function call(element::Element, field_name::ASCIIString, ip, time, ::Type{Val{:G
element(ip, time, Val{:Grad})*element[field_name](time)
end
function call(element::Element, field_name::ASCIIString, ip, time=0.0)
field = element[field_name](time)
function call(element::Element, field_name::ASCIIString, time)
return element[field_name](time)
end
function call(element::Element, field_name::ASCIIString, ip, time)
field = element(field_name, time)
isa(field, DCTI) && return field.data
basis = element(ip, time)
n = length(element)
@@ -141,6 +144,11 @@ function get_connectivity(element::Element)
end
function get_integration_points(element::Element)
# first time initialize default integration points
if length(element.integration_points) == 0
ips = get_integration_points(element.properties)
element.integration_points = [IP(i, w, xi) for (i, (w, xi)) in enumerate(ips)]
end
return element.integration_points
end
src/problems.jl
+1 -1
View File
@@ -176,7 +176,7 @@ function update_assembly!(problem, u, la)
elseif get_formulation_type(problem) == :incremental
info("$(problem.name): incremental formulation, adding increment to solution vector")
assembly.u += u
assembly.la += la
assembly.la = la
elseif get_formulation_type(problem) == :forwarddiff
info("$(problem.name): forwarddiff formulation, adding increment to solution vector")
assembly.u += u