From bd1e64c661455f7ec0bf2fc6bee7b96405fbde0a Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Wed, 12 Sep 2018 17:57:01 +0300 Subject: [PATCH] refactor, normals calculation done somehow different .. --- src/contact2dad.jl | 104 ++++++++++++++++++++++++++++++--------------- test/test_02.jl | 38 +++++++++++++---- 2 files changed, 98 insertions(+), 44 deletions(-) diff --git a/src/contact2dad.jl b/src/contact2dad.jl index fac56ec..94ea7fc 100644 --- a/src/contact2dad.jl +++ b/src/contact2dad.jl @@ -170,11 +170,6 @@ function project_from_slave_to_master_ad( end -""" -Frictionless 2d finite sliding contact with forwarddiff. - -true/false flags: finite_sliding, friction, use_forwarddiff -""" function FEMBase.assemble_elements!(problem::Problem{Contact2DAD}, assembly::Assembly, elements::Vector{Element{Seg2}}, time::Float64) @@ -199,30 +194,69 @@ function FEMBase.assemble_elements!(problem::Problem{Contact2DAD}, assembly::Ass C = zeros(T, 2, nnodes) # 1. update nodal normals for slave elements - normals = similar(u) - tangents = similar(u) - nadj_nodes = Dict{Int64, Int64}() - coeff = props.rotate_normals ? -1.0 : 1.0 - for slave_element in slave_elements - c1, c2 = get_connectivity(slave_element) - x1, x2 = x[c1], x[c2] - tangent = (x2-x1) / norm(x2-x1) - normal = [-tangent[2], tangent[1]] - for c in (c1, c2) - if !haskey(nadj_nodes, c) - normals[c] = zeros(2) - tangents[c] = zeros(2) - nadj_nodes[c] = 0 + normals = Dict{Int64, Vector{T}}() + tangents = Dict{Int64, Vector{T}}() + + # nadj_nodes = Dict{Int64, Int64}() + # coeff = props.rotate_normals ? -1.0 : 1.0 + # for slave_element in slave_elements + # c1, c2 = get_connectivity(slave_element) + # x1, x2 = x[c1], x[c2] + # tangent = (x2-x1) / norm(x2-x1) + # normal = [-tangent[2], tangent[1]] + # for c in (c1, c2) + # if !haskey(nadj_nodes, c) + # normals[c] = zeros(T, 2) + # tangents[c] = zeros(T, 2) + # nadj_nodes[c] = 0 + # end + # normals[c] += coeff*normal + # tangents[c] += coeff*tangent + # nadj_nodes[c] += 1 + # end + # end + # for j in keys(normals) + # normals[j] /= nadj_nodes[j] + # tangents[j] /= nadj_nodes[j] + # end + + #normals = empty(u) + #tangents = empty(u) + Q = [0.0 -1.0; 1.0 0.0] + for element in slave_elements + a, b = conn = get_connectivity(element) + X_el = (X[a], X[b]) + x_el = (X[a]+u[a], X[b]+u[b]) + dN = get_dbasis(element, [0.0], time) + t = sum([kron(dN[:,i], x_el[i]') for i=1:length(x_el)]) + n = Q*t' + n /= norm(n) + t /= norm(t) + for c in conn + if !haskey(normals, c) + normals[c] = vec(n) + else + normals[c] += vec(n) + end + if !haskey(tangents, c) + tangents[c] = vec(t) + else + tangents[c] += vec(t) end - normals[c] += coeff*normal - tangents[c] += coeff*tangent - nadj_nodes[c] += 1 end end - for j in keys(normals) - normals[j] /= nadj_nodes[j] - tangents[j] /= nadj_nodes[j] + for nid in keys(normals) + normals[nid] /= norm(normals[nid]) + tangents[nid] /= norm(tangents[nid]) end + # swap element normals in 2d if they point to inside of body + if props.rotate_normals + for j in keys(normals) + normals[j] = -normals[j] + tangents[j] = -tangents[j] + end + end + update!(slave_elements, "normal", time => normals) update!(slave_elements, "tangent", time => tangents) @@ -236,11 +270,11 @@ function FEMBase.assemble_elements!(problem::Problem{Contact2DAD}, assembly::Ass slave_element_nodes = get_connectivity(slave_element) c1, c2 = get_connectivity(slave_element) X1 = (X[c1], X[c2]) - u1 = ((u[i] for i in slave_element_nodes)...) - x1 = ((Xi+ui for (Xi,ui) in zip(X1,u1))...) - la1 = ((la[i] for i in slave_element_nodes)...) - n1 = ((normals[i] for i in slave_element_nodes)...) - t1 = ((tangents[i] for i in slave_element_nodes)...) + u1 = ((u[i] for i in slave_element_nodes)...,) + x1 = ((Xi+ui for (Xi,ui) in zip(X1,u1))...,) + la1 = ((la[i] for i in slave_element_nodes)...,) + n1 = ((normals[i] for i in slave_element_nodes)...,) + t1 = ((tangents[i] for i in slave_element_nodes)...,) nnodes = size(slave_element, 2) Ae = Matrix(1.0I, nnodes, nnodes) @@ -268,8 +302,8 @@ function FEMBase.assemble_elements!(problem::Problem{Contact2DAD}, assembly::Ass master_element_nodes = get_connectivity(master_element) cm1, cm2 = get_connectivity(master_element) X2 = master_element("geometry", time) - u2 = ((u[i] for i in master_element_nodes)...) - x2 = ((Xi+ui for (Xi,ui) in zip(X2,u2))...) + u2 = ((u[i] for i in master_element_nodes)...,) + x2 = ((Xi+ui for (Xi,ui) in zip(X2,u2))...,) # calculate segmentation: we care only about endpoints xi1a = project_from_master_to_slave_ad(slave_element, field(x1), field(n1), x2[1]) @@ -313,8 +347,8 @@ function FEMBase.assemble_elements!(problem::Problem{Contact2DAD}, assembly::Ass master_element_nodes = get_connectivity(master_element) X2 = master_element("geometry", time) - u2 = ((u[i] for i in master_element_nodes)...) - x2 = ((Xi+ui for (Xi,ui) in zip(X2,u2))...) + u2 = ((u[i] for i in master_element_nodes)...,) + x2 = ((Xi+ui for (Xi,ui) in zip(X2,u2))...,) #x1_midpoint = 1/2*(x1[1]+x1[2]) #x2_midpoint = 1/2*(x2[1]+x2[2]) @@ -411,7 +445,7 @@ function FEMBase.assemble_elements!(problem::Problem{Contact2DAD}, assembly::Ass @info("Summary of nodes") for j in sort(collect(keys(is_active))) n = map(ForwardDiff.value, normals[j]) - info("$j, c=$(condition[j]), s=$(is_active[j]), n=$n") + @info("$j, c=$(condition[j]), s=$(is_active[j]), n=$n") end end diff --git a/test/test_02.jl b/test/test_02.jl index f6d1a48..400a60f 100644 --- a/test/test_02.jl +++ b/test/test_02.jl @@ -41,22 +41,42 @@ function to_dict(u, ndofs, nnodes) return Dict(j => [u[ndofs*(j-1)+i] for i=1:ndofs] for j=1:nnodes) end +function unpack_full(assembly, ndofs) + K = Matrix(assembly.K, ndofs, ndofs) + C1 = Matrix(assembly.C1, ndofs, ndofs) + C2 = Matrix(assembly.C2, ndofs, ndofs) + D = Matrix(assembly.D, ndofs, ndofs) + f = Vector(assembly.f, ndofs) + g = Vector(assembly.g, ndofs) + return K, C1, C2, D, f, g +end + +function calc_err(A, B) + err = norm(A - B) + err_rel = err / max(norm(A), norm(B)) + return err, err_rel +end + ndofs = 20 contact.properties.iteration = 1 -data = i1 -time = 0.0 for (data, time) in zip([i1, i2, i3], [0.0, 1.0, 2.0]) - println("Testing data for iteration / time $time ") empty!(contact.assembly) contact.assembly.u = vec(data["u"]) contact.assembly.la = vec(data["la"]) update!(contact, "displacement", time => to_dict(data["u"], 2, 5)) update!(contact, "lambda", time => to_dict(data["la"], 2, 5)) assemble!(contact, time) - @test isapprox(Matrix(sparse(contact.assembly.K, ndofs, ndofs)), data["K"]) - @test isapprox(Matrix(sparse(contact.assembly.C1, ndofs, ndofs)), data["C1"]) - @test isapprox(Matrix(sparse(contact.assembly.C2, ndofs, ndofs)), data["C2"]) - @test isapprox(Matrix(sparse(contact.assembly.D, ndofs, ndofs)), data["D"]) - @test isapprox(Vector(sparse(contact.assembly.f, ndofs, 1)[:]), data["f"]; atol=1.0e-9) - @test isapprox(Vector(sparse(contact.assembly.g, ndofs, 1)[:]), data["g"]) + + K, C1, C2, D, f, g = unpack_full(contact.assembly, ndofs) + Kerr, Kerr_rel = calc_err(K, data["K"]) + C2err, C2err_rel = calc_err(C2, data["C2"]) + Derr, Derr_rel = calc_err(D, data["D"]) + @debug "Iteration at time $time" Kerr Kerr_rel C2err C2err_rel Derr Derr_rel + + @test isapprox(K, data["K"]) + @test isapprox(C1, data["C1"]) + @test isapprox(C2, data["C2"]) + @test isapprox(D, data["D"]) + @test isapprox(f, data["f"]; atol=1.0e-9) + @test isapprox(g, data["g"]) end