From 3ac6e6924ee6a31ce73e0cef35ef0adf7c223177 Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Mon, 4 Jul 2016 02:52:42 +0300 Subject: [PATCH] frictionless 3d small sliding contact working again --- src/JuliaFEM.jl | 1 + src/problems_contact_2d.jl | 25 ++- src/problems_contact_3d.jl | 283 ++++++++++++++++++++++++++ test/test_contact_2d_small_sliding.jl | 3 +- test/test_contact_3d_small_sliding.jl | 60 ++++++ 5 files changed, 361 insertions(+), 11 deletions(-) create mode 100644 src/problems_contact_3d.jl create mode 100644 test/test_contact_3d_small_sliding.jl diff --git a/src/JuliaFEM.jl b/src/JuliaFEM.jl index c348fdf..f03d5ca 100644 --- a/src/JuliaFEM.jl +++ b/src/JuliaFEM.jl @@ -87,6 +87,7 @@ export calculate_normals, ### Mortar methods, contact mechanics extension ### include("problems_contact.jl") include("problems_contact_2d.jl") +include("problems_contact_3d.jl") export Contact module API diff --git a/src/problems_contact_2d.jl b/src/problems_contact_2d.jl index adc63f9..52b505a 100644 --- a/src/problems_contact_2d.jl +++ b/src/problems_contact_2d.jl @@ -3,7 +3,16 @@ typealias ContactElements2D Union{Seg2} -""" Frictionless 2d small sliding contact without forwarddiff. """ +""" +Frictionless 2d small sliding contact without forwarddiff. + +problem +time +dimension +finite_sliding +friction +use_forwarddiff +""" function assemble!(problem::Problem{Contact}, time::Float64, ::Type{Val{1}}, ::Type{Val{false}}, ::Type{Val{false}}, ::Type{Val{false}}; debug=false) @@ -23,11 +32,11 @@ function assemble!(problem::Problem{Contact}, time::Float64, for slave_element in slave_elements nsl = length(slave_element) - X1 = slave_element["geometry"](time) - u1 = slave_element["displacement"](time) - la1 = slave_element["reaction force"](time) - n1 = slave_element["normal"](time) - t1 = slave_element["tangent"](time) + X1 = slave_element("geometry", time) + u1 = slave_element("displacement", time) + la1 = slave_element("reaction force", time) + n1 = slave_element("normal", time) + t1 = slave_element("tangent", time) x1 = X1 + u1 Q1_ = [n1[1] t1[1]] Q2_ = [n1[2] t1[2]] @@ -38,7 +47,7 @@ function assemble!(problem::Problem{Contact}, time::Float64, if "element area" in props.store_fields element_area = 0.0 - for ip in get_integration_points(slave_element, 3) + for ip in get_integration_points(slave_element) detJ = slave_element(ip, time, Val{:detJ}) w = ip.weight*detJ element_area += w @@ -243,6 +252,4 @@ function assemble!(problem::Problem{Contact}, time::Float64, problem.assembly.D = D problem.assembly.g = g - return - end diff --git a/src/problems_contact_3d.jl b/src/problems_contact_3d.jl new file mode 100644 index 0000000..2598c62 --- /dev/null +++ b/src/problems_contact_3d.jl @@ -0,0 +1,283 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md + +typealias ContactElements3D Union{Tri3, Tri6, Quad4, Quad8, Quad9} + +function create_orthogonal_basis(n) + I = eye(3) + k = indmax([norm(cross(n,I[:,k])) for k in 1:3]) + t1 = cross(n, I[:,k])/norm(cross(n, I[:,k])) + t2 = cross(n, t1) + return t1, t2 +end + +""" +Frictionless 2d small sliding contact. + +problem +time +dimension +finite_sliding +friction +use_forwarddiff +""" +function assemble!(problem::Problem{Contact}, time::Float64, + ::Type{Val{2}}, ::Type{Val{false}}, + ::Type{Val{false}}, ::Type{Val{false}}; debug=true) + + props = problem.properties + field_dim = get_unknown_field_dimension(problem) + field_name = get_parent_field_name(problem) + slave_elements = get_slave_elements(problem) + + # 1. calculate nodal normals and tangents for slave element nodes j ∈ S + normals = calculate_normals(slave_elements, time, Val{2}; + rotate_normals=props.rotate_normals) + update!(slave_elements, "normal", normals) + + # 2. loop all slave elements + for (slave_num, slave_element) in enumerate(slave_elements) + + nsl = length(slave_element) + X1 = slave_element("geometry", time) + u1 = slave_element("displacement", time) + la = slave_element("reaction force", time) + n1 = slave_element("normal", time) + t11, t21 = create_orthogonal_basis(n1[1]) + t12, t22 = create_orthogonal_basis(n1[2]) + t13, t23 = create_orthogonal_basis(n1[3]) + Q1_ = [n1[1] t11 t21] + Q2_ = [n1[2] t12 t22] + Q3_ = [n1[3] t13 t23] + Z = zeros(3, 3) + Q3 = [Q1_ Z Z; Z Q2_ Z; Z Z Q3_] + contact_area = 0.0 + contact_error = 0.0 + + element_area = 0.0 + for ip in get_integration_points(slave_element) + detJ = slave_element(ip, time, Val{:detJ}) + w = ip.weight*detJ + element_area += w + end + if "element area" in props.store_fields + update!(slave_element, "element area", time => element_area) + end + + if slave_num == 1 + info("First slave element area = $element_area") + info("NT basis of first slave element") + dump(Q3) + end + + # project slave nodes to auxiliary plane (x0, Q) + #xi = get_reference_element_midpoint(slave_element) + xi = [1/3, 1/3] + N = vec(get_basis(slave_element, xi, time)) + x0 = N*X1 + n0 = N*n1 + S = Vector[project_vertex_to_auxiliary_plane(p, x0, n0) for p in X1] + + # 3. loop all master elements + for master_element in slave_element("master elements", time) + + nm = length(master_element) + X2 = master_element("geometry", time) + 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 + norm(mean(X1) - X2[3]) / norm(X1[2] - X1[1]) < props.distval || continue + =# + + # 3.1 project master nodes to auxiliary plane and create polygon clipping + M = Vector[project_vertex_to_auxiliary_plane(p, x0, n0) for p in X2] + P = get_polygon_clip(S, M, n0) + length(P) < 3 && continue # no clipping or shared edge (no volume) + check_orientation!(P, n0) + C0 = calculate_centroid(P) + + De = zeros(nsl, nsl) + Me = zeros(nsl, nm) + ge = zeros(field_dim*nsl) + + # 4. loop integration cells + for cell in get_cells(P, C0) + virtual_element = Element(Tri3) + update!(virtual_element, "geometry", cell) + + # 5. loop integration point of integration cell + for ip in get_integration_points(virtual_element, 3) + + # project gauss point from auxiliary plane to master and slave element + x_gauss = virtual_element("geometry", ip, time) + if isnan(x_gauss[1]) + info("is nan") + info("x_gauss = $x_gauss") + info("cell = $cell") + info("C0 = $C0") + info("P = $P") + info("S = $S") + info("M = $M") + info("n0 = $n0") + error("nan, unable to continue") + end + xi_s, alpha = project_vertex_to_surface(x_gauss, x0, n0, slave_element, X1, time) + xi_m, alpha = project_vertex_to_surface(x_gauss, x0, n0, master_element, X2, time) + + detJ = virtual_element(ip, time, Val{:detJ}) + w = ip.weight*detJ + # add contributions + N1 = vec(get_basis(slave_element, xi_s, time)) + N2 = vec(get_basis(master_element, xi_m, time)) + De += w*N1*N1' + Me += w*N1*N2' + + x_s = N1*(X1+u1) + x_m = N2*(X2+u2) + ge += w*vec((x_m-x_s)*N1') + contact_area += w + n_s = N1*n1 + contact_error += 1/2*w*dot(n_s, x_s-x_m)^2 + end # integration points done + + end # integration cells done + + # 6. add contribution to contact virtual work + sdofs = get_gdofs(problem, slave_element) + mdofs = get_gdofs(problem, master_element) + nsldofs = length(sdofs) + nmdofs = length(mdofs) + D3 = zeros(nsldofs, nsldofs) + M3 = zeros(nmdofs, nmdofs) + for i=1:field_dim + D3[i:field_dim:end, i:field_dim:end] += De + M3[i:field_dim:end, i:field_dim:end] += Me + end + + add!(problem.assembly.C1, sdofs, sdofs, D3) + add!(problem.assembly.C1, sdofs, mdofs, -M3) + add!(problem.assembly.C2, sdofs, sdofs, Q3'*D3) + add!(problem.assembly.C2, sdofs, mdofs, -Q3'*M3) + add!(problem.assembly.g, sdofs, Q3'*ge) + + end # master elements done + + if "contact area" in props.store_fields + update!(slave_element, "contact area", time => contact_area) + end + + end # slave elements done, contact virtual work ready + + S = sort(collect(keys(normals))) # slave element nodes + weighted_gap = Dict{Int64, Vector{Float64}}() + contact_pressure = Dict{Int64, Vector{Float64}}() + complementarity_condition = Dict{Int64, Vector{Float64}}() + is_active = Dict{Int64, Int}() + is_inactive = Dict{Int64, Int}() + is_slip = Dict{Int64, Int}() + is_stick = Dict{Int64, Int}() + + g = full(problem.assembly.g) + la = problem.assembly.la + + # active / inactive node detection + for j in S + dofs = [3*(j-1)+1, 3*(j-1)+2, 3*(j-1)+3] + weighted_gap[j] = g[dofs] + if length(la) != 0 + normal = normals[j] + tangent1, tangent2 = create_orthogonal_basis(normal) + p = dot(normal, la[dofs]) + t1 = dot(tangent1, la[dofs]) + t2 = dot(tangent2, la[dofs]) + contact_pressure[j] = [p, t1, t2] + else + contact_pressure[j] = [0.0, 0.0, 0.0] + end + complementarity_condition[j] = contact_pressure[j] - weighted_gap[j] + if complementarity_condition[j][1] < 0 + is_inactive[j] = 1 + is_active[j] = 0 + is_slip[j] = 0 + is_stick[j] = 0 + else + is_inactive[j] = 0 + is_active[j] = 1 + is_slip[j] = 1 + is_stick[j] = 0 + end + end + + if "weighted gap" in props.store_fields + update!(slave_elements, "weighted gap", time => weighted_gap) + end + if "contact pressure" in props.store_fields + update!(slave_elements, "contact pressure", time => contact_pressure) + end + if "complementarity condition" in props.store_fields + update!(slave_elements, "complementarity condition", time => complementarity_condition) + end + if "active nodes" in props.store_fields + update!(slave_elements, "active nodes", time => is_active) + end + if "inactive nodes" in props.store_fields + update!(slave_elements, "inactive nodes", time => is_inactive) + end + if "stick nodes" in props.store_fields + update!(slave_elements, "stick nodes", time => is_stick) + end + if "slip nodes" in props.store_fields + update!(slave_elements, "slip nodes", time => is_slip) + end + + info("# | active | inactive | stick | slip | gap | pres | comp") + for j in S + str1 = "$j | $(is_active[j]) | $(is_inactive[j]) | $(is_stick[j]) | $(is_slip[j]) | " + str2 = "$(round(weighted_gap[j], 3)) | $(round(contact_pressure[j], 3)) | $(round(complementarity_condition[j], 3))" + info(str1 * str2) + end + + # 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 = [3*(j-1)+1, 3*(j-1)+2, 3*(j-1)+3] + tdofs = dofs[[2,3]] + if (is_active[j] == 1) && (is_slip[j] == 1) + info("$j is in active/slip, removing tangential constraints $tdofs") + C2[tdofs,:] = 0.0 + g[tdofs] = 0.0 + normal = normals[j] + tangent1, tangent2 = create_orthogonal_basis(normal) + D[tdofs[1], dofs] = tangent1 + D[tdofs[2], dofs] = tangent2 + end + end + + # remove inactive nodes from assembly + for j in S + dofs = [3*(j-1)+1, 3*(j-1)+2, 3*(j-1)+3] + if is_inactive[j] == 1 + info("$j is inactive, removing dofs $dofs") + C1[dofs,:] = 0.0 + C2[dofs,:] = 0.0 + D[dofs,:] = 0.0 + g[dofs,:] = 0.0 + end + end + + problem.assembly.C1 = C1 + problem.assembly.C2 = C2 + problem.assembly.D = D + problem.assembly.g = g + +end + diff --git a/test/test_contact_2d_small_sliding.jl b/test/test_contact_2d_small_sliding.jl index 62406e7..4d8ae33 100644 --- a/test/test_contact_2d_small_sliding.jl +++ b/test/test_contact_2d_small_sliding.jl @@ -45,7 +45,6 @@ function get_model(::Type{Val{Symbol("curved 2d contact small sliding")}}) interface_master_elements = create_elements(mesh, "UPPER_BOTTOM") update!(interface_slave_elements, "master elements", interface_master_elements) interface.elements = [interface_master_elements; interface_slave_elements] - info("type of list is ", typeof(first(interface_slave_elements)("master elements", 0.0))) solver = Solver(Nonlinear) push!(solver, upper, lower, bc_upper, bc_lower, interface) @@ -70,7 +69,7 @@ end function get_model(::Type{Val{Symbol("hertz contact, full 2d model")}}) # from fenet d3613 advanced finite element contact benchmarks # a = 6.21 mm, pmax = 3585 MPa - # this is a very dense mesh and for that reason pmax is not very + # this is a very sparse mesh and for that reason pmax is not very # (only 6 elements in -20 .. 20 mm contact zone, 3 elements in contact # instead integrate pressure in normal and tangential direction mesh = get_mesh("hertz contact, full 2d model") diff --git a/test/test_contact_3d_small_sliding.jl b/test/test_contact_3d_small_sliding.jl new file mode 100644 index 0000000..54997fd --- /dev/null +++ b/test/test_contact_3d_small_sliding.jl @@ -0,0 +1,60 @@ +# 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.Test + +@testset "3d upper side curved contact" begin + + # TODO: accurate solution is not known, verify using another fem software + # however results look very meaningful and probably this is right. + + meshfile = Pkg.dir("JuliaFEM") * "/test/testdata/block_3d_curved.med" + mesh = aster_read_mesh(meshfile) + + upper = Problem(Elasticity, "upper", 3) + upper.elements = create_elements(mesh, "UPPER") + update!(upper, "youngs modulus", 96.0) + update!(upper, "poissons ratio", 1/3) + + lower = Problem(Elasticity, "lower", 3) + lower.elements = create_elements(mesh, "LOWER") + update!(lower, "youngs modulus", 96.0) + update!(lower, "poissons ratio", 1/3) + + bc_upper = Problem(Dirichlet, "upper boundary", 3, "displacement") + bc_upper.elements = create_elements(mesh, "UPPER_TOP") + update!(bc_upper, "displacement 1", 0.0) + update!(bc_upper, "displacement 2", 0.0) + update!(bc_upper, "displacement 3", -0.1) + + bc_lower = Problem(Dirichlet, "lower boundary", 3, "displacement") + bc_lower.elements = create_elements(mesh, "LOWER_BOTTOM") + update!(bc_lower, "displacement 1", 0.0) + update!(bc_lower, "displacement 2", 0.0) + update!(bc_lower, "displacement 3", 0.0) + + contact = Problem(Contact, "contact between upper and lower block", 3, "displacement") + contact_slave_elements = create_elements(mesh, "LOWER_TOP") + contact_master_elements = create_elements(mesh, "UPPER_BOTTOM") + update!(contact_slave_elements, "master elements", contact_master_elements) + contact.elements = [contact_master_elements; contact_slave_elements] + + solver = Solver(Nonlinear) + push!(solver, upper, lower, bc_upper, bc_lower, contact) + + call(solver) + for element in get_slave_elements(contact) + normal = element("normal", [1/3, 1/3], solver.time) + @test isapprox(normal, [0.0, 0.0, 1.0]) + pres = dot(normal, element("reaction force", [1/3, 1/3], solver.time)) + info("pressure = $pres") + #info(element("displacement", [1/3, 1/3], solver.time)) + end + normu = norm(contact.assembly.u) + info("displacement field norm = $normu") + @test isapprox(normu, 0.7417557629004985) + +end