mirror of
https://github.com/JuliaFEM/JuliaFEM.jl.git
synced 2026-09-27 04:10:54 +00:00
Forwarddiff elasticity working again.
This commit is contained in:
Jukka Aho
+110
-81
@@ -6,10 +6,11 @@ type Elasticity <: FieldProblem
|
||||
# these are found from problem.properties for type Problem{Elasticity}
|
||||
formulation :: Symbol
|
||||
finite_strain :: Bool
|
||||
use_forwarddiff :: Bool
|
||||
end
|
||||
function Elasticity()
|
||||
# formulations: plane_stress, plane_strain, continuum
|
||||
return Elasticity(:continuum, true)
|
||||
return Elasticity(:continuum, true, false)
|
||||
end
|
||||
|
||||
# in case of experimenting new things;
|
||||
@@ -31,7 +32,9 @@ end
|
||||
function assemble!(assembly::Assembly, problem::Problem{Elasticity}, element::Element, time::Real)
|
||||
props = problem.properties
|
||||
gdofs = get_gdofs(problem, element)
|
||||
if props.formulation == :continuum
|
||||
if props.use_forwarddiff
|
||||
Kt, f = assemble(problem, element, time, Val{:forwarddiff})
|
||||
elseif props.formulation == :continuum
|
||||
Kt, f = assemble(problem, element, time, Val{:continuum})
|
||||
elseif (props.formulation == :plane_stress) || (props.formulation == :plane_strain)
|
||||
Kt, f = assemble(problem, element, time, Val{:plane})
|
||||
@@ -306,6 +309,111 @@ function assemble{El<:Union{Tri3, Tri6, Quad4}}(problem::Problem{Elasticity}, el
|
||||
return Kt, f
|
||||
end
|
||||
|
||||
""" Elasticity equations using ForwardDiff
|
||||
|
||||
Formulation
|
||||
-----------
|
||||
|
||||
Field equation is:
|
||||
∂u/∂t = ∇⋅f - b
|
||||
|
||||
Weak form is: find u∈U such that ∀v in V
|
||||
|
||||
δW := ∫ρ₀∂²u/∂t²⋅δu dV₀ + ∫S:δE dV₀ - ∫b₀⋅δu dV₀ - ∫t₀⋅δu dA₀ = 0
|
||||
|
||||
where
|
||||
|
||||
ρ₀ = density
|
||||
b₀ = displacement load
|
||||
t₀ = displacement traction
|
||||
|
||||
References
|
||||
----------
|
||||
|
||||
https://en.wikipedia.org/wiki/Linear_elasticity
|
||||
https://en.wikipedia.org/wiki/Finite_strain_theory
|
||||
https://en.wikipedia.org/wiki/Stress_measures
|
||||
https://en.wikipedia.org/wiki/Mooney%E2%80%93Rivlin_solid
|
||||
https://en.wikipedia.org/wiki/Strain_energy_density_function
|
||||
https://en.wikipedia.org/wiki/Plane_stress
|
||||
https://en.wikipedia.org/wiki/Hooke's_law
|
||||
|
||||
"""
|
||||
function assemble(problem::Problem{Elasticity}, element::Element, time::Real, ::Type{Val{:forwarddiff}})
|
||||
|
||||
dim = get_unknown_field_dimension(problem)
|
||||
nnodes = size(element, 2)
|
||||
|
||||
function get_residual_vector(u::Vector)
|
||||
u = reshape(u, dim, nnodes)
|
||||
u = Field([u[:,i] for i=1:nnodes])
|
||||
r = zeros(dim, nnodes)
|
||||
|
||||
for ip in get_integration_points(element)
|
||||
|
||||
JT = transpose(get_jacobian(element, ip, time))
|
||||
n, m = size(JT)
|
||||
if n == m
|
||||
w = ip.weight*det(JT)
|
||||
elseif m == 1
|
||||
w = ip.weight*norm(JT)
|
||||
elseif m == 2
|
||||
w = ip.weight*norm(cross(JT[:,1], JT[:,2]))
|
||||
else
|
||||
error("jacobian $JT")
|
||||
end
|
||||
|
||||
# calculate internal forces
|
||||
if haskey(element, "youngs modulus") && haskey(element, "poissons ratio")
|
||||
grad = element(ip, time, Val{:grad})
|
||||
gradu = grad*u
|
||||
|
||||
# kinematics
|
||||
F = I + gradu
|
||||
E = 1/2*(F'*F - I)
|
||||
|
||||
# material
|
||||
young = element("youngs modulus", ip, time)
|
||||
poisson = element("poissons ratio", ip, time)
|
||||
mu = young/(2*(1+poisson))
|
||||
lambda = young*poisson/((1+poisson)*(1-2*poisson))
|
||||
if problem.properties.formulation == :plane_stress
|
||||
lambda = 2*lambda*mu/(lambda + 2*mu) # <- correction for plane stress
|
||||
end
|
||||
|
||||
# stress
|
||||
S = lambda*trace(E)*I + 2*mu*E
|
||||
|
||||
r += w*F*S*grad
|
||||
end
|
||||
|
||||
# calculate external forces - volume load
|
||||
if haskey(element, "displacement load")
|
||||
basis = element(ip, time)
|
||||
b = element("displacement load", ip, time)
|
||||
r -= w*b*basis
|
||||
end
|
||||
|
||||
# external forces - surface traction force
|
||||
if haskey(element, "displacement traction force")
|
||||
basis = element(ip, time)
|
||||
T = element("displacement traction force", ip, time)
|
||||
r -= w*T*basis
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
return vec(r)
|
||||
|
||||
end
|
||||
|
||||
field = element("displacement", time)
|
||||
Kt, allresults = ForwardDiff.jacobian(get_residual_vector, vec(field),
|
||||
AllResults, cache=autodiffcache)
|
||||
f = -ForwardDiff.value(allresults)
|
||||
return Kt, f
|
||||
end
|
||||
|
||||
|
||||
###############################
|
||||
# Plastic material #
|
||||
@@ -376,84 +484,5 @@ function get_unknown_field_type{P<:ElasticityProblem}(::Type{P})
|
||||
end
|
||||
|
||||
|
||||
""" Elasticity equations.
|
||||
|
||||
Formulation
|
||||
-----------
|
||||
|
||||
Field equation is:
|
||||
∂u/∂t = ∇⋅f - b
|
||||
|
||||
Weak form is: find u∈U such that ∀v in V
|
||||
|
||||
δW := ∫ρ₀∂²u/∂t²⋅δu dV₀ + ∫S:δE dV₀ - ∫b₀⋅δu dV₀ - ∫t₀⋅δu dA₀ = 0
|
||||
|
||||
where
|
||||
|
||||
ρ₀ = density
|
||||
b₀ = displacement load
|
||||
t₀ = displacement traction
|
||||
|
||||
References
|
||||
----------
|
||||
|
||||
https://en.wikipedia.org/wiki/Linear_elasticity
|
||||
https://en.wikipedia.org/wiki/Finite_strain_theory
|
||||
https://en.wikipedia.org/wiki/Stress_measures
|
||||
https://en.wikipedia.org/wiki/Mooney%E2%80%93Rivlin_solid
|
||||
https://en.wikipedia.org/wiki/Strain_energy_density_function
|
||||
https://en.wikipedia.org/wiki/Plane_stress
|
||||
https://en.wikipedia.org/wiki/Hooke's_law
|
||||
|
||||
"""
|
||||
function get_residual_vector{P<:ElasticityProblem}(problem::Problem{P}, element::Element, ip::IntegrationPoint, time::Number; variation=nothing)
|
||||
r = zeros(Float64, problem.dim, length(element))
|
||||
|
||||
J = get_jacobian(element, ip, time)
|
||||
|
||||
# internal forces
|
||||
if haskey(element, "youngs modulus") && haskey(element, "poissons ratio")
|
||||
u = element("displacement", time, variation)
|
||||
grad = element(ip, time, Val{:grad})
|
||||
gradu = grad*u
|
||||
|
||||
# deformation gradient
|
||||
F = I + gradu
|
||||
|
||||
# material
|
||||
young = element("youngs modulus", ip, time)
|
||||
poisson = element("poissons ratio", ip, time)
|
||||
mu = young/(2*(1+poisson))
|
||||
lambda = young*poisson/((1+poisson)*(1-2*poisson))
|
||||
if problem.properties.formulation == :plane_stress
|
||||
lambda = 2*lambda*mu/(lambda + 2*mu) # <- correction for 2d problems
|
||||
end
|
||||
|
||||
# strain
|
||||
E = 1/2*(F'*F - I)
|
||||
# stress
|
||||
S = lambda*trace(E)*I + 2*mu*E
|
||||
|
||||
r += F*S*grad*det(J)
|
||||
end
|
||||
|
||||
# external forces - volume load
|
||||
if haskey(element, "displacement load")
|
||||
basis = element(ip, time)
|
||||
b = element("displacement load", ip, time)
|
||||
r -= b*basis*det(J)
|
||||
end
|
||||
|
||||
# external forces - surface traction force
|
||||
if haskey(element, "displacement traction force")
|
||||
basis = element(ip, time)
|
||||
T = element("displacement traction force", ip, time)
|
||||
JT = transpose(J)
|
||||
s = size(JT, 2) == 1 ? JT : cross(JT[:,1], JT[:,2])
|
||||
r -= T*basis*norm(s)
|
||||
end
|
||||
|
||||
return vec(r)
|
||||
end
|
||||
|
||||
=#
|
||||
|
||||
+2
-3
@@ -279,7 +279,7 @@ function Base.call(field::CCTV, time::Number)
|
||||
return field.data(time)
|
||||
end
|
||||
|
||||
### Interpolation
|
||||
### Interpolation
|
||||
|
||||
""" Interpolate time-invariant field in time direction. """
|
||||
function Base.call(field::DVTI, time::Float64)
|
||||
@@ -361,11 +361,10 @@ function Base.call(basis::CVTI, xi::Vector, time::Number)
|
||||
call(basis, xi)
|
||||
end
|
||||
|
||||
function Base.(:*)(grad::Matrix{Float64}, field::DVTI)
|
||||
function Base.(:*)(grad::Matrix, field::DVTI)
|
||||
return sum([kron(grad[:,i], field[i]') for i=1:length(field)])'
|
||||
end
|
||||
|
||||
### FIELDSET ###
|
||||
|
||||
typealias FieldSet Dict{ASCIIString, Field}
|
||||
|
||||
|
||||
+68
-32
@@ -1,46 +1,82 @@
|
||||
# This file is a part of JuliaFEM.
|
||||
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
|
||||
|
||||
module ElasticityTests
|
||||
|
||||
using JuliaFEM.Test
|
||||
using JuliaFEM
|
||||
using JuliaFEM.Core: Seg2, Quad4, Hex8,
|
||||
ElasticityProblem, PlaneStressElasticityProblem,
|
||||
solve!, get_connectivity, DirichletProblem
|
||||
using JuliaFEM.Core: Node, Seg2, Quad4, Elasticity, Dirichlet, Problem, Solver, update!
|
||||
using JuliaFEM.Core: assemble
|
||||
|
||||
|
||||
function test_elasticity_volume_load()
|
||||
@testset "test forwarddiff version + volume load." begin
|
||||
nodes = Dict{Int64, Node}(
|
||||
1 => [0.0, 0.0],
|
||||
2 => [10.0, 0.0],
|
||||
3 => [10.0, 1.0],
|
||||
4 => [0.0, 1.0])
|
||||
# constant volume load on nodes
|
||||
load = Dict(
|
||||
1 => [0.0, -10.0],
|
||||
2 => [0.0, -10.0],
|
||||
3 => [0.0, -10.0],
|
||||
4 => [0.0, -10.0])
|
||||
young = Dict(1 => 500.0, 2 => 500.0, 3 => 500.0, 4 => 500.0)
|
||||
poisson = Dict(1 => 0.3, 2 => 0.3, 3 => 0.3, 4 => 0.3)
|
||||
element = Quad4([1, 2, 3, 4])
|
||||
element["geometry"] = Vector[[0.0, 0.0], [10.0, 0.0], [10.0, 1.0], [0.0, 1.0]]
|
||||
element["youngs modulus"] = 500.0
|
||||
element["poissons ratio"] = 0.3
|
||||
element["displacement load"] = Vector[[0.0, -10.0], [0.0, -10.0], [0.0, -10.0], [0.0, -10.0]]
|
||||
element["displacement"] = (0.0 => Vector{Float64}[[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]])
|
||||
problem = PlaneStressElasticityProblem()
|
||||
push!(problem, element)
|
||||
update!(element, "geometry", nodes)
|
||||
update!(element, "youngs modulus", young)
|
||||
update!(element, "poissons ratio", poisson)
|
||||
update!(element, "displacement load", load)
|
||||
boundary = Seg2([1, 4])
|
||||
update!(boundary, "geometry", nodes)
|
||||
update!(boundary, "displacement 1", 0.0)
|
||||
update!(boundary, "displacement 2", 0.0)
|
||||
|
||||
free_dofs = [3, 4, 5, 6]
|
||||
solve!(problem, free_dofs, 0.0; max_iterations=10)
|
||||
body = Problem(Elasticity, "beam", 2)
|
||||
body.properties.formulation = :plane_stress
|
||||
body.properties.use_forwarddiff = true
|
||||
push!(body, element)
|
||||
bc = Problem(Dirichlet, "fixed left side", 2, "displacement")
|
||||
#bc.properties.formulation = :incremental
|
||||
push!(bc, boundary)
|
||||
|
||||
solver = Solver()
|
||||
push!(solver, body, bc)
|
||||
call(solver)
|
||||
disp = element("displacement", [1.0, 1.0], 0.0)
|
||||
# function get_previous_ip(element::Element, current_ip::IntegrationPoint)
|
||||
# end
|
||||
# ipdata = element("integration points", time) => IntegrationPoint[ip1, ip2, ..., ipN]
|
||||
# for some_ip in ipdata
|
||||
# if isapprox(some_ip.xi, ip.xi)
|
||||
# info("found")
|
||||
# last_value = some_ip("material parameter", time)
|
||||
# break
|
||||
# end
|
||||
# end
|
||||
#ip1 = last(element["integration points"])[1]
|
||||
#ip2 = last(element["integration points"])[2]
|
||||
# strain = ip1("gl strain")
|
||||
info("displacement at tip: $disp")
|
||||
#info("strain in first ip: $strain. ip coord = $(ip1.xi) and weight = $(ip1.weight)")
|
||||
# verified using Code Aster, verification/2015-10-22-plane-stress/cplan_grot_gdep_volume_force.resu
|
||||
@test isapprox(disp[2], -8.77303119819776)
|
||||
end
|
||||
#test_elasticity_volume_load()
|
||||
|
||||
@testset "test that stiffness matrix is same" begin
|
||||
nodes = Dict{Int64, Node}(
|
||||
1 => [0.0, 0.0],
|
||||
2 => [10.0, 0.0],
|
||||
3 => [10.0, 1.0],
|
||||
4 => [0.0, 1.0])
|
||||
displacement = Dict(
|
||||
1 => [0.1, 0.2],
|
||||
2 => [0.3, 0.4],
|
||||
3 => [0.5, 0.6],
|
||||
4 => [0.7, 0.8])
|
||||
displacement = Dict(
|
||||
1 => [0.0, 0.0],
|
||||
2 => [0.0, 0.0],
|
||||
3 => [0.0, 0.0],
|
||||
4 => [0.0, 0.0])
|
||||
load = Dict(
|
||||
1 => [0.0, -10.0],
|
||||
2 => [0.0, -10.0],
|
||||
3 => [0.0, -10.0],
|
||||
4 => [0.0, -10.0])
|
||||
element = Quad4([1, 2, 3, 4])
|
||||
update!(element, "geometry", nodes)
|
||||
update!(element, "displacement", displacement)
|
||||
update!(element, "youngs modulus", 288.0)
|
||||
update!(element, "poissons ratio", 1/3)
|
||||
update!(element, "displacement load", load)
|
||||
body = Problem(Elasticity, "beam", 2)
|
||||
body.properties.formulation = :plane_stress
|
||||
K1, f1 = assemble(body, element, 0.0, Val{:forwarddiff})
|
||||
K2, f2 = assemble(body, element, 0.0, Val{:plane})
|
||||
@test isapprox(K1, K2)
|
||||
@test isapprox(f1, f2)
|
||||
end
|
||||
|
||||
Reference in New Issue
Block a user