diff --git a/src/elasticity.jl b/src/elasticity.jl index 035388b..b564aef 100644 --- a/src/elasticity.jl +++ b/src/elasticity.jl @@ -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 =# diff --git a/src/fields.jl b/src/fields.jl index 38857e0..61681bd 100644 --- a/src/fields.jl +++ b/src/fields.jl @@ -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} - diff --git a/test/test_elasticity.jl b/test/test_elasticity.jl index 8249df2..f3c0e06 100644 --- a/test/test_elasticity.jl +++ b/test/test_elasticity.jl @@ -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