From 2173034c2f1bfb5b0589651d4afba84878156f2b Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Mon, 12 Sep 2016 03:39:33 +0300 Subject: [PATCH] minor changes to 2d contact formulation --- src/elements.jl | 20 +++---- src/problems_contact.jl | 3 +- src/problems_contact_2d.jl | 87 ++++++++++++++++++++-------- src/problems_mortar_2d.jl | 11 ++-- test/test_elasticity_plane_strain.jl | 45 ++++++++++++++ 5 files changed, 127 insertions(+), 39 deletions(-) create mode 100644 test/test_elasticity_plane_strain.jl diff --git a/src/elements.jl b/src/elements.jl index 98d43ab..577b6f4 100644 --- a/src/elements.jl +++ b/src/elements.jl @@ -158,19 +158,19 @@ As a result element now have time invariant (variable) vector field "geometry" w """ function update!(element::Element, field_name, data::Dict) - element[field_name] = Field(data) - #element[field_name] = [data[i] for i in get_connectivity(element)] + #element[field_name] = Field(data) + element[field_name] = [data[i] for i in get_connectivity(element)] end function update!{K,V}(element::Element, field_name, data::Pair{Float64, Dict{K, V}}) - #time, field_data = data - #element_data = V[field_data[i] for i in get_connectivity(element)] - #update!(element, field_name, time => element_data) - if haskey(element, field_name) - update!(element[field_name], data) - else - element[field_name] = Field(data) - end + time, field_data = data + element_data = V[field_data[i] for i in get_connectivity(element)] + update!(element, field_name, time => element_data) + #if haskey(element, field_name) + # update!(element[field_name], data) + #else + # element[field_name] = Field(data) + #end end function update!(element::Element, field_name::AbstractString, datas::Union{Real, Vector, Pair{Float64, Union{Float64, Real, Vector{Any}}}}...) diff --git a/src/problems_contact.jl b/src/problems_contact.jl index 4adbe40..f15abef 100644 --- a/src/problems_contact.jl +++ b/src/problems_contact.jl @@ -19,6 +19,7 @@ type Contact <: BoundaryProblem use_forwarddiff :: Bool minimum_active_set_size :: Int distval :: Float64 + remove_from_set :: Bool # allow removal of non-potential contact pairs store_fields :: Vector{AbstractString} end @@ -26,7 +27,7 @@ function Contact() default_fields = ["element area", "contact area", "weighted gap", "contact pressure", "active nodes", "inactive nodes", "stick nodes", "slip nodes", "complementarity condition", "contact error"] - return Contact(-1, false, false, false, true, false, 0, 5.0, default_fields) + return Contact(-1, false, false, false, true, false, 0, 5.0, false, default_fields) end function get_unknown_field_name(problem::Problem{Contact}) diff --git a/src/problems_contact_2d.jl b/src/problems_contact_2d.jl index 20fc8a0..e117b6e 100644 --- a/src/problems_contact_2d.jl +++ b/src/problems_contact_2d.jl @@ -21,6 +21,8 @@ function assemble!(problem::Problem{Contact}, time::Float64, update!(slave_elements, "normal", time => normals) update!(slave_elements, "tangent", time => tangents) + Rn = 0.0 + # 2. loop all slave elements for slave_element in slave_elements @@ -56,8 +58,13 @@ function assemble!(problem::Problem{Contact}, time::Float64, u2 = master_element("displacement", time) x2 = X2 + u2 - norm(mean(X1) - X2[1]) / norm(X1[2] - X1[1]) < props.distval || continue - norm(mean(X1) - X2[2]) / norm(X1[2] - X1[1]) < props.distval || continue + if norm(mean(X1) - X2[1]) / norm(X1[2] - X1[1]) > props.distval + continue + end + + if norm(mean(X1) - X2[2]) / norm(X1[2] - X1[1]) > props.distval + continue + end # 3.1 calculate segmentation xi1a = project_from_master_to_slave(slave_element, X2[1], time) @@ -89,7 +96,11 @@ function assemble!(problem::Problem{Contact}, time::Float64, # local mortar matrices fill!(De, 0.0) fill!(Me, 0.0) - ge = zeros(field_dim*nsl) + Ne = zeros(nsl, 2*nsl) + Te = zeros(nsl, 2*nsl) + He = zeros(nsl, 2*nsl) + ce = zeros(nsl) + ge = zeros(nsl) for ip in get_integration_points(slave_element, 3) detJ = slave_element(ip, time, Val{:detJ}) w = ip.weight*detJ*l @@ -97,40 +108,56 @@ function assemble!(problem::Problem{Contact}, time::Float64, xi_s = dot([1/2*(1-xi); 1/2*(1+xi)], xi1) N1 = vec(get_basis(slave_element, xi_s, time)) Phi = Ae*N1 + # project gauss point from slave element to master element in direction n_s X_s = N1*X1 # coordinate in gauss point n_s = N1*n1 # normal direction in gauss point + t_s = N1*t1 # tangent condition in gauss point + n_s /= norm(n_s) + t_s /= norm(t_s) xi_m = project_from_slave_to_master(master_element, X_s, n_s, time) N2 = vec(get_basis(master_element, xi_m, time)) X_m = N2*X2 - De += w*Phi*N1' - Me += w*Phi*N2' + u_s = N1*u1 u_m = N2*u2 x_s = X_s + u_s x_m = X_m + u_m la_s = Phi*la1 - ge += w*vec((x_m-x_s)*Phi') + + # virtual work + De += w*Phi*N1' + Me += w*Phi*N2' + + # contact constraints + Ne += w*reshape(kron(N1, n_s, Phi), 2, 4) + Te += w*reshape(kron(N2, n_s, Phi), 2, 4) + He += w*reshape(kron(N1, t_s, Phi), 2, 4) + ge += w*Phi*dot(n_s, x_m-x_s) + ce += w*N1*dot(n_s, -la_s) + Rn += w*dot(n_s, -la_s) + contact_area += w contact_error += 1/2*w*dot(n_s, x_s-x_m)^2 end - # add contribution to contact virtual work sdofs = get_gdofs(problem, slave_element) mdofs = get_gdofs(problem, master_element) - nsldofs = length(sdofs) - nmdofs = length(mdofs) - D2 = zeros(nsldofs, nsldofs) - M2 = zeros(nmdofs, nmdofs) + + # add contribution to contact virtual work for i=1:field_dim - D2[i:field_dim:end, i:field_dim:end] += De - M2[i:field_dim:end, i:field_dim:end] += Me + lsdofs = sdofs[i:field_dim:end] + lmdofs = mdofs[i:field_dim:end] + add!(problem.assembly.C1, lsdofs, lsdofs, De) + add!(problem.assembly.C1, lsdofs, lmdofs, -Me) end - add!(problem.assembly.C1, sdofs, sdofs, D2) - add!(problem.assembly.C1, sdofs, mdofs, -M2) - add!(problem.assembly.C2, sdofs, sdofs, Q2'*D2) - add!(problem.assembly.C2, sdofs, mdofs, -Q2'*M2) - add!(problem.assembly.g, sdofs, Q2'*ge) + + # add contribution to contact constraints + add!(problem.assembly.C2, sdofs[1:field_dim:end], sdofs, Ne) + add!(problem.assembly.C2, sdofs[1:field_dim:end], mdofs, -Te) + add!(problem.assembly.D, sdofs[2:field_dim:end], sdofs, He) + add!(problem.assembly.g, sdofs[1:field_dim:end], ge) + add!(problem.assembly.c, sdofs[1:field_dim:end], ce) end # master elements done @@ -153,13 +180,22 @@ function assemble!(problem::Problem{Contact}, time::Float64, is_slip = Dict{Int64, Int}() is_stick = Dict{Int64, Int}() - g = full(problem.assembly.g) la = problem.assembly.la + ndofs = length(la) + info("contact ndofs: $ndofs") + info("Rn = $Rn") + + C1 = sparse(problem.assembly.C1) + C2 = sparse(problem.assembly.C2, ndofs, ndofs) + D = sparse(problem.assembly.D, ndofs, ndofs) + g = full(problem.assembly.g, ndofs, 1) + c = full(problem.assembly.c, ndofs, 1) # active / inactive node detection for j in S dofs = [2*(j-1)+1, 2*(j-1)+2] weighted_gap[j] = g[dofs] + if length(la) != 0 p = dot(normals[j], la[dofs]) t = dot(tangents[j], la[dofs]) @@ -167,6 +203,8 @@ function assemble!(problem::Problem{Contact}, time::Float64, else contact_pressure[j] = [0.0, 0.0] end + +# contact_pressure[j] = c[dofs] complementarity_condition[j] = contact_pressure[j] - weighted_gap[j] if complementarity_condition[j][1] < 0 is_inactive[j] = 1 @@ -178,6 +216,11 @@ function assemble!(problem::Problem{Contact}, time::Float64, is_active[j] = 1 is_slip[j] = 1 is_stick[j] = 0 + _c1 = complementarity_condition[j][1] + _c2 = c[dofs] + _c3 = contact_pressure[j][1] + _c4 = g[dofs] + info("active $j: c1 = $_c1, c2 = $_c2, c3 = $_c3, c4 = $_c4") end end @@ -212,12 +255,9 @@ function assemble!(problem::Problem{Contact}, time::Float64, # solve variational inequality - C1 = sparse(problem.assembly.C1) - ndofs = size(C1, 1) - C2 = sparse(problem.assembly.C2) - D = spzeros(ndofs, ndofs) # constitutive modelling in tangent direction, frictionless contact + #= for j in S dofs = [2*(j-1)+1, 2*(j-1)+2] if (is_active[j] == 1) && (is_slip[j] == 1) @@ -227,6 +267,7 @@ function assemble!(problem::Problem{Contact}, time::Float64, D[dofs[2], dofs] = tangents[j] end end + =# # remove inactive nodes from assembly for j in S diff --git a/src/problems_mortar_2d.jl b/src/problems_mortar_2d.jl index 7934bd4..7ddf890 100644 --- a/src/problems_mortar_2d.jl +++ b/src/problems_mortar_2d.jl @@ -68,9 +68,10 @@ function calculate_normals(elements, time, ::Type{Val{1}}; rotate_normals=false) tangents = Dict{Int64, Vector{Float64}}() for element in elements conn = get_connectivity(element) - X1 = element("geometry", time) - dN = get_dbasis(element, [0.0], time) - tangent = vec(sum([kron(dN[:,i], X1[i]') for i=1:length(X1)])) + #X1 = element("geometry", time) + #dN = get_dbasis(element, [0.0], time) + #tangent = vec(sum([kron(dN[:,i], X1[i]') for i=1:length(X1)])) + tangent = vec(element([0.0], time, Val{:Jacobian})) for nid in conn if haskey(tangents, nid) tangents[nid] += tangent @@ -116,8 +117,8 @@ function assemble!(problem::Problem{Mortar}, time::Float64, ::Type{Val{1}}, ::Ty # 1. calculate nodal normals and tangents for slave element nodes j ∈ S normals, tangents = calculate_normals(slave_elements, time, Val{1}; rotate_normals=props.rotate_normals) - update!(slave_elements, "normal", normals) - update!(slave_elements, "tangent", tangents) + update!(slave_elements, "normal", time => normals) + update!(slave_elements, "tangent", time => tangents) # 2. loop all slave elements for slave_element in slave_elements diff --git a/test/test_elasticity_plane_strain.jl b/test/test_elasticity_plane_strain.jl new file mode 100644 index 0000000..33ed666 --- /dev/null +++ b/test/test_elasticity_plane_strain.jl @@ -0,0 +1,45 @@ +# 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.Postprocess +using JuliaFEM.Testing + +@testset "test 2d linear elasticity with surface + volume load" begin + meshfile = "/geometry/2d_block/BLOCK_1elem.med" + mesh = aster_read_mesh(Pkg.dir("JuliaFEM")*meshfile) + + # field problem + block = Problem(Elasticity, "BLOCK", 2) + block.properties.formulation = :plane_strain + block.properties.finite_strain = false + block.properties.geometric_stiffness = false + block.elements = create_elements(mesh, "BLOCK") + update!(block.elements, "youngs modulus", 288.0) + update!(block.elements, "poissons ratio", 1/3) + + # traction + traction = Problem(Elasticity, "TRACTION", 2) + traction.properties.formulation = :plane_strain + traction.properties.finite_strain = false + traction.properties.geometric_stiffness = false + traction.elements = create_elements(mesh, "TOP") + update!(traction, "displacement traction force 2", 288.0*9/8) + + # boundary conditions + bc_sym_23 = Problem(Dirichlet, "symmetry bc 23", 2, "displacement") + bc_sym_23.elements = create_elements(mesh, "LEFT") + update!(bc_sym_23, "displacement 1", 0.0) + bc_sym_13 = Problem(Dirichlet, "symmetry bc 13", 2, "displacement") + bc_sym_13.elements = create_elements(mesh, "BOTTOM") + update!(bc_sym_13, "displacement 2", 0.0) + + solver = LinearSolver(block, traction, bc_sym_23, bc_sym_13) + solver() + + info("u = ", block.assembly.u) + info("λ = ", block.assembly.la) + +end +