diff --git a/src/mortar_2d_autodiff.jl b/src/mortar_2d_autodiff.jl index b4d2cbb..7b696f3 100644 --- a/src/mortar_2d_autodiff.jl +++ b/src/mortar_2d_autodiff.jl @@ -82,416 +82,9 @@ function project_from_slave_to_master{E<:MortarElements2D}( end error("find projection from slave to master: did not converge, last val: $xi2 and $dxi2") - -end - - -function assemble!{E<:MortarElements2D}(assembly::Assembly, problem::Problem{Mortar}, slave_element::Element{E}, time::Real, - ::Type{Val{:forwarddiff_old}}) - haskey(slave_element, "master elements") || return - props = problem.properties - field_dim = get_unknown_field_dimension(problem) - field_name = get_parent_field_name(problem) - slave_dofs = get_gdofs(slave_element, field_dim) - nnodes = size(slave_element, 2) - X1 = slave_element("geometry", time) - #u1 = slave_element("displacement", time) - #x1 = X1 + u1 - slave_element_nodes = get_connectivity(slave_element) - adjacent_elements = find_elements(get_elements(problem), slave_element_nodes) - adjacent_nodes = get_nodes(adjacent_elements) # including also nodes from adjacent elements - Q = [0.0 -1.0; 1.0 0.0] - - X = spzeros(10000, 1) - for element in get_elements(problem) - conn = get_connectivity(element) - geom = element("geometry", time) - for (c, g) in zip(conn, geom) - dofs = [field_dim*(c-1)+1, field_dim*(c-1)+2] - X[dofs] = g - end - end - - # here x does not mean deformed configuration - x = [problem.assembly.u; problem.assembly.la] - if length(x) == 0 - info("mortar_2d_autodiff: length(x) == 0") - # resize solution vectors according to initial configuration of this problem - X = vec(full(sparse(findnz(X)...))) - x = zeros(length(X)*2) - else - # resize initial configuration to match real dimension - I, J, V = findnz(X) - X = vec(full(sparse(I, J, V, length(problem.assembly.u), 1))) - end - ndofs = round(Int, length(x)/2) - # at the end we should have -# info("mortar_2d_autodiff: size of x = $(size(x))") -# info("mortar_2d_autodiff: size of X = $(size(X))") -# info("mortar_2d_autodiff: ndofs = $ndofs") - - """ Calculate normal vector for slave element nodes in current configuration. """ - function calculate_normals(u::Matrix) - normals = zeros(u) - # 1. update nodal normals - for element in adjacent_elements - conn = get_connectivity(element) - gdofs = get_gdofs(element, field_dim) - X_el = element("geometry", time) - u_el = Field(Vector[u[:, i] for i in conn]) - x_el = X_el + u_el - for ip in get_integration_points(element, Val{3}) - dN = get_dbasis(element, ip) - N = element(ip, time) - t = sum([kron(dN[:,i], x_el[i]') for i=1:length(x_el)]) - normals[:, conn] += ip.weight*Q*t'*N - end - end - slave_normals = Field(Vector[normals[:,i]/norm(normals[:,i]) for i in slave_element_nodes]) - return slave_normals - end - - function calculate_mortar_projection(u::Matrix, n1::DVTI) - B = SparseMatrixCOO{Real}([], [], []) - - u1 = Field([u[:,i] for i in slave_element_nodes]) - x1 = X1 + u1 - - for master_element in slave_element["master elements"] - X2 = master_element("geometry", time) - master_element_nodes = get_connectivity(master_element) - u2 = Field([u[:,i] for i in master_element_nodes]) - x2 = X2 + u2 - #info("master element coordinate 1 = $(ForwardDiff.get_value(x2[1]))") - #info("master element coordinate 2 = $(ForwardDiff.get_value(x2[2]))") - - # calculate segmentation: we care only about endpoints - # note: these are quadratic/cubic functions, analytical solution possible - xi1a = project_from_master_to_slave(slave_element, x1, n1, x2[1]) - xi1b = project_from_master_to_slave(slave_element, x1, n1, x2[end]) - xi1 = clamp([xi1a; xi1b], -1.0, 1.0) - l = 1/2*abs(xi1[2]-xi1[1]) - isapprox(l, 0.0) && continue # no contribution - - # integrate slave side - D = zeros(nnodes, nnodes) - Me = zeros(nnodes, nnodes) - for ip in get_integration_points(slave_element, Val{5}) - dN = get_dbasis(slave_element, ip) - # jacobian of slave element in deformed state - j = sum([kron(dN[:,i], x1[i]') for i=1:length(x1)]) - w = ip.weight*norm(j)*l - xi_s = dot([1/2*(1-ip.xi); 1/2*(1+ip.xi)], xi1) - N = get_basis(slave_element, xi_s) - D += w*diagm(vec(N)) - Me += w*N'*N - end - Ae = D*inv(Me) - - # integrate master side - M = zeros(nnodes, nnodes) - for ip in get_integration_points(slave_element, Val{5}) - dN = get_dbasis(slave_element, ip) - # jacobian of slave element in deformed state - j = sum([kron(dN[:,i], x1[i]') for i=1:length(x1)]) - w = ip.weight*norm(j)*l - xi_g = dot([1/2*(1-ip.xi); 1/2*(1+ip.xi)], xi1) - N1 = get_basis(slave_element, xi_g) - x_g = vec(N1)*x1 - n_g = vec(N1)*n1 - #info("slave gauss point coordinates $(ForwardDiff.get_value(x_g))") - #info("slave gauss point normal direction $(ForwardDiff.get_value(n_g))") - xi_m = project_from_slave_to_master(master_element, x_g, n_g, x2) - N2 = get_basis(master_element, xi_m) - M += w*kron(Ae*N1', N2) - end - - slave_dofs = get_gdofs(slave_element, field_dim) - master_dofs = get_gdofs(master_element, field_dim) - for i=1:field_dim - add!(B, slave_dofs[i:field_dim:end], slave_dofs[i:field_dim:end], D) - add!(B, slave_dofs[i:field_dim:end], master_dofs[i:field_dim:end], -M) - end - - end - - return B - end - - function calculate_contact_rhs(x::Vector) - ndofs = round(Int, length(x)/2) - u = x[1:ndofs] - la = x[ndofs+1:end] -# info("calculate_contact_rhs: size of u = $(size(u))") -# info("calculate_contact_rhs: size of X = $(size(X))") -# info("calculate_contact_rhs: size of la = $(size(la))") -# info("calculate_contact_rhs: ndofs = $ndofs") - u2 = reshape(u, field_dim, round(Int, length(u)/field_dim)) - - normals = calculate_normals(u2) - B = calculate_mortar_projection(u2, normals) - B = sparse(B, ndofs, ndofs) - fc = B' * la # contact force residual for r = fint + fc - fext = 0 - - N = SparseMatrixCOO{Real}([], [], []) - T = SparseMatrixCOO{Real}([], [], []) - for (i, j) in enumerate(slave_element_nodes) - dofs = [2*(j-1)+1, 2*(j-1)+2] - add!(N, dofs, [dofs[1]], reshape(normals[i], 2, 1)) - add!(T, dofs, [dofs[2]], reshape(Q'*normals[i], 2, 1)) - end - N = sparse(N, ndofs, ndofs) - T = sparse(T, ndofs, ndofs) - gn = -N*B*(X+u) - gt = -T*B*(X+u) -# gn = -N*B*u - lan = N*la - lat = T*la - cn = 1.0e3 - C = lan - max(0, lan - cn*gn) + lat - # C = lan - gn + lat - gt <-- ihan viturallensa - # C = gn+gt <-- not working - # C = B*(X+u) - # C = B*u <- pitää kiinni, "tie". - # C = N*B*u + T*la <- palikat menee väärään suuntaan - # C = -N*B*(X+u) + T*la <- toimii suht hyvin mut kääntyy väärään suuntaan (t-suunnassa) - # C = -N*B*(X+u) - T*la <- sama - # C = -N*B*(X+u) - T*B*la <- sama - # C = N*la + T*la - max(0, N*la + N*B*(X+u)) - cond = lan[1:field_dim:end] - cn*gn[1:field_dim:end] - all_nodes = slave_element_nodes - inactive_nodes = find(cond .<= 0) - active_nodes = find(cond .> 0) - inactive_nodes = setdiff(all_nodes, inactive_nodes) - active_nodes = setdiff(all_nodes, active_nodes) - info("S = $all_nodes, I = $inactive_nodes, A = $active_nodes") - info("lambda = $(ForwardDiff.get_value(lan[slave_dofs]))") - info("gn = $(ForwardDiff.get_value(gn[slave_dofs]))") - #for j in active_nodes - # dofs = [field_dim*(j-1)+1, field_dim*(j-1)+2] - # C[dofs] = 0 - #end - -# C = -(N+T)*B*(X+u) -# C = -B*(X+u) - - return [fc; C] - - end - - A, allresults = ForwardDiff.jacobian(calculate_contact_rhs, x, ForwardDiff.AllResults) - b = -ForwardDiff.value(allresults) - A = sparse(A) - b = sparse(b) - K = A[1:ndofs,1:ndofs] - C1 = A[1:ndofs,ndofs+1:end]' - C2 = A[ndofs+1:end,1:ndofs] - D = A[ndofs+1:end,ndofs+1:end] - f = b[1:ndofs] - g = b[ndofs+1:end] - C2[2:field_dim:end] = 0 - g[2:field_dim:end] = 0 - -# C2[2:field_dim:end] = 0 -# g[2:field_dim:end] = 0 -# inactives = find(g[1:field_dim:end] .<= 0) -# actives = find(g[1:field_dim:end] .> 0) -# inactives = setdiff(slave_element_nodes, inactives) -# actives = setdiff(slave_element_nodes, actives) -# info("all nodes = $slave_element_nodes, inactives = $inactives, actives = $actives") -# info("g = $g") - -# for j in inactives - # dofs = [field_dim*(j-1)+1, field_dim*(j-1)+2] - # K[dofs,:] = 0 - # C1[dofs,:] = 0 - # C2[dofs,:] = 0 - # D[dofs,:] = 0 - # g[dofs,:] = 0 - #end - #for j in actives - # dofs = [field_dim*(j-1)+1, field_dim*(j-1)+2] - # C2[dofs[1],:] = 0 - #end - add!(assembly.K, K) - add!(assembly.C1, C1) - add!(assembly.C2, C2) - add!(assembly.D, D) - add!(assembly.f, f) - add!(assembly.g, g) -end - - -function assemble!{E<:MortarElements2D}(assembly::Assembly, - problem::Problem{Mortar}, slave_element::Element{E}, - time::Real, ::Type{Val{:forwarddiff_old2}}) - haskey(slave_element, "master elements") || return - props = problem.properties - field_dim = get_unknown_field_dimension(problem) - field_name = get_parent_field_name(problem) - - function calculate_interface(x::Vector) - - ndofs = round(Int, length(x)/2) - nnodes = round(Int, ndofs/field_dim) - u = reshape(x[1:ndofs], field_dim, nnodes) - la = reshape(x[ndofs+1:end], field_dim, nnodes) - fc = zeros(u) - C = zeros(la) - - slave_element_nodes = get_connectivity(slave_element) - X1 = slave_element("geometry", time) - u1 = Field(Vector[u[:,i] for i in slave_element_nodes]) - la1 = Field(Vector[la[:,i] for i in slave_element_nodes]) - x1 = X1 + u1 - - # 1. update nodal normals for this element. average nodes from adjacent elements - adjacent_elements = find_elements(get_elements(problem), slave_element_nodes) - adjacent_nodes = get_nodes(adjacent_elements) # including also nodes from adjacent elements - Q = [0.0 -1.0; 1.0 0.0] - normals = zeros(u) - for element in adjacent_elements - conn = get_connectivity(element) - gdofs = get_gdofs(element, field_dim) - X_el = element("geometry", time) - u_el = Field(Vector[u[:, i] for i in conn]) - x_el = X_el + u_el - for ip in get_integration_points(element, Val{3}) - dN = get_dbasis(element, ip) - N = element(ip, time) - t = sum([kron(dN[:,i], x_el[i]') for i=1:length(x_el)]) - normals[:, conn] += ip.weight*Q*t'*N - end - end - # --> slave side normals in deformed state - n1 = Field(Vector[normals[:,i]/norm(normals[:,i]) for i in slave_element_nodes]) - - for master_element in slave_element["master elements"] - - master_element_nodes = get_connectivity(master_element) - X2 = master_element("geometry", time) - u2 = Field(Vector[u[:,i] for i in master_element_nodes]) - x2 = X2 + u2 - - # calculate segmentation: we care only about endpoints - # note: these are quadratic/cubic functions, analytical solution possible - xi1a = project_from_master_to_slave(slave_element, x1, n1, x2[1]) - xi1b = project_from_master_to_slave(slave_element, x1, n1, x2[end]) - xi1 = clamp([xi1a; xi1b], -1.0, 1.0) - l = 1/2*abs(xi1[2]-xi1[1]) - isapprox(l, 0.0) && continue # no contribution in this master element - - nnodes = size(slave_element, 2) - De = zeros(nnodes, nnodes) - Me = zeros(nnodes, nnodes) - for ip in get_integration_points(slave_element, Val{5}) - # jacobian of slave element in deformed state - dN = get_dbasis(slave_element, ip) - j = sum([kron(dN[:,i], x1[i]') for i=1:length(x1)]) - w = ip.weight*norm(j)*l - xi_s = dot([1/2*(1-ip.xi); 1/2*(1+ip.xi)], xi1) - N1 = get_basis(slave_element, xi_s) - De += w*diagm(vec(N1)) - Me += w*N1'*N1 - end - Ae = De*inv(Me) - - slave_dofs = get_gdofs(slave_element, field_dim) - master_dofs = get_gdofs(master_element, field_dim) - - for ip in get_integration_points(slave_element, Val{5}) - # jacobian of slave element in deformed state - dN = get_dbasis(slave_element, ip) - j = sum([kron(dN[:,i], x1[i]') for i=1:length(x1)]) - w = ip.weight*norm(j)*l - - # project gauss point from slave element to master element - xi_s = dot([1/2*(1-ip.xi); 1/2*(1+ip.xi)], xi1) - N1 = vec(get_basis(slave_element, xi_s)) - x_s = N1*x1 # coordinate in gauss point - n_s = N1*n1 # normal direction in gauss point - t_s = Q'*n_s # tangent direction in gauss point - R_s = [n_s t_s] - xi_m = project_from_slave_to_master(master_element, x_s, n_s, x2) - N2 = vec(get_basis(master_element, xi_m)) - x_m = N2*x2 - Phi = Ae*N1 - - u_s = N1*u1 - u_m = N2*u2 - - la_s = Phi*la1 # traction force in gauss point - #lan = dot(n_s, la_s) # normal component - #lat = dot(t_s, la_s) # tangential component - la_nt = R_s*la_s - g = x_s-x_m - gn = props.gap_sign*dot(n_s, g) # normal gap - #gu = dot(n_s, u_s - u_m) # normal displacement gap - #gt = dot(t_s, x_s - x_m) # tangential gap - - fc[:,slave_element_nodes] += w*la_s*N1' - fc[:,master_element_nodes] -= w*la_s*N2' - C[1,slave_element_nodes] += w*gn*Phi' - #C[2,slave_element_nodes] += w*la_nt[2]*Phi' - #C[2,slave_element_nodes] += w*la_nt[2,:]*Phi' - #C[2,slave_element_nodes] += w*dot(t_s, u_s - u_m)*Phi' # <-- for tie - - #R += w*R_s - - end - - end # master elements done - - for (i, j) in enumerate(slave_element_nodes) - n = n1[i] - t = Q'*n - R = [n t] - la_nt = R*la[:,j] - C[2,j] += la_nt[2] - end - - return vec([fc C]) - - end - - # x doesn't mean deformed configuration here - x = [problem.assembly.u; problem.assembly.la] - ndofs = round(Int, length(x)/2) - #if ndofs == 0 - # info("INITIALIZING THINGS") - # problem.assembly.u = zeros(16) - # problem.assembly.la = zeros(16) - # x = [problem.assembly.u; problem.assembly.la] - # ndofs = round(Int, length(x)/2) - #end - - A, allresults = ForwardDiff.jacobian(calculate_interface, x, ForwardDiff.AllResults) - b = -ForwardDiff.value(allresults) - #b = -calculate_interface(x) - #info("PE = $(ForwardDiff.value(allresults))") - A = sparse(A) - b = sparse(b) - SparseMatrix.droptol!(A, 1.0e-12) - SparseMatrix.droptol!(b, 1.0e-12) - #println(A) - K = A[1:ndofs,1:ndofs] - C1 = transpose(A[1:ndofs,ndofs+1:end]) - C2 = A[ndofs+1:end,1:ndofs] - D = A[ndofs+1:end,ndofs+1:end] - f = b[1:ndofs] - g = b[ndofs+1:end] - add!(assembly.K, K) - add!(assembly.C1, C1) - add!(assembly.C2, C2) - add!(assembly.D, D) - add!(assembly.f, f) - add!(assembly.g, g) - - return end - function assemble!{E<:MortarElements2D}(assembly::Assembly, problem::Problem{Mortar}, slave_element::Element{E}, time::Real, ::Type{Val{:forwarddiff}}) @@ -508,7 +101,7 @@ function assemble!{E<:MortarElements2D}(assembly::Assembly, la = reshape(x[ndofs+1:end], field_dim, nnodes) fc = zeros(u) C = zeros(la) - + slave_element_nodes = get_connectivity(slave_element) X1 = slave_element("geometry", time) u1 = Field(Vector[u[:,i] for i in slave_element_nodes]) @@ -601,7 +194,7 @@ function assemble!{E<:MortarElements2D}(assembly::Assembly, fc[:,slave_element_nodes] += w*la_s*N1' fc[:,master_element_nodes] -= w*la_s*N2' #C[1,slave_element_nodes] += w*gn*Phi' - + lan_tot += w*la_nt[1]*Phi gap_tot += w*gn*Phi @@ -609,15 +202,6 @@ function assemble!{E<:MortarElements2D}(assembly::Assembly, end # master elements done -# ncf = lan_tot - max(0, lan_tot - gap_tot) -# info("pressure in nodes: $(ForwardDiff.get_value(lan_tot))") -# info("weighted gap in nodes: $(ForwardDiff.get_value(gap_tot))") -# info("ncf: $(ForwardDiff.get_value(ncf))") -# cond = +lan_tot + gap_tot -# cond = +lan_tot - gap_tot # singular -# cond = -lan_tot + gap_tot -# cond = -lan_tot - gap_tot - for (i, j) in enumerate(slave_element_nodes) n = n1[i] t = Q'*n @@ -627,9 +211,8 @@ function assemble!{E<:MortarElements2D}(assembly::Assembly, # if -lan_tot[i] + gap_tot[i] < 0 if -la_nt[1] + gap_tot[i] < 0 - #if j in [31, 32, 33, 34, 35, 36] info("set node $j active") - C[1,j] -= gap_tot[i] + C[1,j] += gap_tot[i] # C[1,j] += la_nt[1] - max(0, la_nt[1] - gap_tot[i]) C[2,j] += la_nt[2] else @@ -671,4 +254,3 @@ function assemble!{E<:MortarElements2D}(assembly::Assembly, return end - diff --git a/src/problems.jl b/src/problems.jl index 3097436..991d5b1 100644 --- a/src/problems.jl +++ b/src/problems.jl @@ -177,7 +177,7 @@ function update_assembly!(problem, u, la) elseif get_formulation_type(problem) == :forwarddiff info("$(problem.name): forwarddiff formulation, adding increment to solution vector") assembly.u += u - assembly.la += la + assembly.la = la else info("$(problem.name): unknown formulation type, don't know what to do with results") error("serious failure with problem formulation: $(get_formulation_type(problem))") diff --git a/test/test_mortar_autodiff.jl b/test/test_mortar_autodiff.jl new file mode 100644 index 0000000..9be705d --- /dev/null +++ b/test/test_mortar_autodiff.jl @@ -0,0 +1,48 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md + +using JuliaFEM.Test + +using JuliaFEM.Core: Node, Seg2, update!, Problem, Mortar, assemble! + +function get_testproblems(u, la) + nodes = Dict{Int64, Node}( + 1 => [0.0, 0.0], + 2 => [2.0, 0.0], + 3 => [1.0, 2.0], + 4 => [0.0, 2.0], + 5 => [2.5, 0.0], + 6 => [4.5, 0.0], + 7 => [4.5, 1.0], + 8 => [2.5, 1.0]) + displacement = Dict{Int64, Vector{Float64}}() + reaction_force = Dict{Int64, Vector{Float64}}() + for i=1:8 + displacement[i] = u[:,i] + reaction_force[i] = la[:,i] + end + bc5 = Seg2([3, 2]) + bc6 = Seg2([8, 5]) + update!([bc5, bc6], "geometry", nodes) + update!([bc5, bc6], "displacement", displacement) + update!([bc5, bc6], "reaction force", reaction_force) + bc5["master elements"] = [bc6] + contact1 = Problem(Mortar, "contact between bodies", 2, "displacement") + contact2 = Problem(Mortar, "contact between bodies", 2, "displacement") + contact2.properties.formulation = :forwarddiff + contact2.assembly.u = vec(u) + contact2.assembly.la = vec(la) + push!(contact1, bc5, bc6) + push!(contact2, bc5, bc6) + return contact1, contact2 +end + +@testset "test linearization of contact force in undeformed state" begin + u = zeros(2, 8) + la = zeros(2, 8) + contact1, contact2 = get_testproblems(u, la) + assemble!(contact1, 0.0) + assemble!(contact2, 0.0) + @test isapprox(full(contact1.assembly.C1), full(contact2.assembly.C1)) + @test isapprox(full(contact1.assembly.K), full(contact2.assembly.K)) +end