From ccbd2d7224cdc4639658fa7584b97a141c510c80 Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Fri, 27 May 2016 17:45:12 +0300 Subject: [PATCH] fixed bug related to strain component ordering and internal force vector --- src/JuliaFEM.jl | 3 +- src/dirichlet.jl | 3 +- src/elasticity.jl | 71 ++++++++++++++++++++++++++++------------------- src/elements.jl | 24 ++++++++++------ src/problems.jl | 2 +- 5 files changed, 63 insertions(+), 40 deletions(-) diff --git a/src/JuliaFEM.jl b/src/JuliaFEM.jl index a085e25..9683a93 100644 --- a/src/JuliaFEM.jl +++ b/src/JuliaFEM.jl @@ -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 diff --git a/src/dirichlet.jl b/src/dirichlet.jl index 6ffb74b..bfdf8c5 100644 --- a/src/dirichlet.jl +++ b/src/dirichlet.jl @@ -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}) diff --git a/src/elasticity.jl b/src/elasticity.jl index e354190..7c84ca3 100644 --- a/src/elasticity.jl +++ b/src/elasticity.jl @@ -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 diff --git a/src/elements.jl b/src/elements.jl index 1142470..9b9a85b 100644 --- a/src/elements.jl +++ b/src/elements.jl @@ -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 diff --git a/src/problems.jl b/src/problems.jl index b7688e8..8cfd39a 100644 --- a/src/problems.jl +++ b/src/problems.jl @@ -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