From 31cd7163bb883e770d0ff540e05d6a3e3dbc79bc Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Mon, 4 Jul 2016 00:08:22 +0300 Subject: [PATCH] removed obsolete code --- src/JuliaFEM.jl | 1 + src/problems_contact.jl | 246 ------------ src/problems_contact_2d.jl | 248 +++++++++++++ src/problems_mortar_3d_autodiff.jl | 513 -------------------------- src/problems_mortar_forwarddiff_pe.jl | 446 ---------------------- 5 files changed, 249 insertions(+), 1205 deletions(-) create mode 100644 src/problems_contact_2d.jl delete mode 100644 src/problems_mortar_3d_autodiff.jl delete mode 100644 src/problems_mortar_forwarddiff_pe.jl diff --git a/src/JuliaFEM.jl b/src/JuliaFEM.jl index 2f102db..c348fdf 100644 --- a/src/JuliaFEM.jl +++ b/src/JuliaFEM.jl @@ -86,6 +86,7 @@ export calculate_normals, ### Mortar methods, contact mechanics extension ### include("problems_contact.jl") +include("problems_contact_2d.jl") export Contact module API diff --git a/src/problems_contact.jl b/src/problems_contact.jl index b6fd476..bdd3c4f 100644 --- a/src/problems_contact.jl +++ b/src/problems_contact.jl @@ -37,8 +37,6 @@ function get_formulation_type(problem::Problem{Contact}) return :incremental end -typealias ContactElements2D Union{Seg2} - function assemble!(problem::Problem{Contact}, time::Real) if problem.properties.dimension == -1 problem.properties.dimension = dim = size(first(problem.elements), 1) @@ -52,247 +50,3 @@ function assemble!(problem::Problem{Contact}, time::Real) assemble!(problem, time, dimension, finite_sliding, friction, use_forwarddiff) end -""" Frictionless 2d small sliding contact without forwarddiff. """ -function assemble!(problem::Problem{Contact}, time::Float64, - ::Type{Val{1}}, ::Type{Val{false}}, - ::Type{Val{false}}, ::Type{Val{false}}; debug=false) - - 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, tangents = calculate_normals(slave_elements, time, Val{1}; - rotate_normals=props.rotate_normals) - update!(slave_elements, "normal", time => normals) - update!(slave_elements, "tangent", time => tangents) - - # 2. loop all slave elements - 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 = X1 + u1 - Q1_ = [n1[1] t1[1]] - Q2_ = [n1[2] t1[2]] - Z = zeros(2, 2) - Q2 = [Q1_ Z; Z Q2_] - contact_area = 0.0 - contact_error = 0.0 - - if "element area" in props.store_fields - element_area = 0.0 - for ip in get_integration_points(slave_element, 3) - detJ = slave_element(ip, time, Val{:detJ}) - w = ip.weight*detJ - element_area += w - end - update!(slave_element, "element area", time => element_area) - end - - # 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 - - # 3.1 calculate segmentation - xi1a = project_from_master_to_slave(slave_element, X2[1], time) - xi1b = project_from_master_to_slave(slave_element, X2[2], time) - 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 - - # 3.2. bi-orthogonal basis - De = zeros(nsl, nsl) - Me = zeros(nsl, nsl) - Ae = zeros(nsl, nsl) - if props.dual_basis - for ip in get_integration_points(slave_element, 3) - detJ = slave_element(ip, time, Val{:detJ}) - w = ip.weight*detJ*l - xi = ip.coords[1] - xi_s = dot([1/2*(1-xi); 1/2*(1+xi)], xi1) - N1 = vec(get_basis(slave_element, xi_s, time)) - De += w*diagm(N1) - Me += w*N1*N1' - end - Ae = De*inv(Me) - else - Ae = eye(nsl) - end - - # 3.3. loop integration points of one integration segment and calculate - # local mortar matrices - fill!(De, 0.0) - fill!(Me, 0.0) - ge = zeros(field_dim*nsl) - for ip in get_integration_points(slave_element, 3) - detJ = slave_element(ip, time, Val{:detJ}) - w = ip.weight*detJ*l - xi = ip.coords[1] - 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 - 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') - 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) - 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 - 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) - - end # master elements done - - if "contact area" in props.store_fields - update!(slave_element, "contact area", time => contact_area) - end - - if "contact error" in props.store_fields - update!(slave_element, "contact error", time => contact_error) - 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 = [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]) - contact_pressure[j] = [p, t] - else - contact_pressure[j] = [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 = [2*(j-1)+1, 2*(j-1)+2] - if (is_active[j] == 1) && (is_slip[j] == 1) - info("$j is in active/slip, removing tangential constraint $(dofs[2])") - C2[dofs[2],:] = 0.0 - g[dofs[2]] = 0.0 - D[dofs[2], dofs] = tangents[j] - end - end - - # remove inactive nodes from assembly - for j in S - dofs = [2*(j-1)+1, 2*(j-1)+2] - 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 - - return - -end - diff --git a/src/problems_contact_2d.jl b/src/problems_contact_2d.jl new file mode 100644 index 0000000..adc63f9 --- /dev/null +++ b/src/problems_contact_2d.jl @@ -0,0 +1,248 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md + +typealias ContactElements2D Union{Seg2} + +""" Frictionless 2d small sliding contact without forwarddiff. """ +function assemble!(problem::Problem{Contact}, time::Float64, + ::Type{Val{1}}, ::Type{Val{false}}, + ::Type{Val{false}}, ::Type{Val{false}}; debug=false) + + 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, tangents = calculate_normals(slave_elements, time, Val{1}; + rotate_normals=props.rotate_normals) + update!(slave_elements, "normal", time => normals) + update!(slave_elements, "tangent", time => tangents) + + # 2. loop all slave elements + 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 = X1 + u1 + Q1_ = [n1[1] t1[1]] + Q2_ = [n1[2] t1[2]] + Z = zeros(2, 2) + Q2 = [Q1_ Z; Z Q2_] + contact_area = 0.0 + contact_error = 0.0 + + if "element area" in props.store_fields + element_area = 0.0 + for ip in get_integration_points(slave_element, 3) + detJ = slave_element(ip, time, Val{:detJ}) + w = ip.weight*detJ + element_area += w + end + update!(slave_element, "element area", time => element_area) + end + + # 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 + + # 3.1 calculate segmentation + xi1a = project_from_master_to_slave(slave_element, X2[1], time) + xi1b = project_from_master_to_slave(slave_element, X2[2], time) + 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 + + # 3.2. bi-orthogonal basis + De = zeros(nsl, nsl) + Me = zeros(nsl, nsl) + Ae = zeros(nsl, nsl) + if props.dual_basis + for ip in get_integration_points(slave_element, 3) + detJ = slave_element(ip, time, Val{:detJ}) + w = ip.weight*detJ*l + xi = ip.coords[1] + xi_s = dot([1/2*(1-xi); 1/2*(1+xi)], xi1) + N1 = vec(get_basis(slave_element, xi_s, time)) + De += w*diagm(N1) + Me += w*N1*N1' + end + Ae = De*inv(Me) + else + Ae = eye(nsl) + end + + # 3.3. loop integration points of one integration segment and calculate + # local mortar matrices + fill!(De, 0.0) + fill!(Me, 0.0) + ge = zeros(field_dim*nsl) + for ip in get_integration_points(slave_element, 3) + detJ = slave_element(ip, time, Val{:detJ}) + w = ip.weight*detJ*l + xi = ip.coords[1] + 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 + 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') + 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) + 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 + 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) + + end # master elements done + + if "contact area" in props.store_fields + update!(slave_element, "contact area", time => contact_area) + end + + if "contact error" in props.store_fields + update!(slave_element, "contact error", time => contact_error) + 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 = [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]) + contact_pressure[j] = [p, t] + else + contact_pressure[j] = [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 = [2*(j-1)+1, 2*(j-1)+2] + if (is_active[j] == 1) && (is_slip[j] == 1) + info("$j is in active/slip, removing tangential constraint $(dofs[2])") + C2[dofs[2],:] = 0.0 + g[dofs[2]] = 0.0 + D[dofs[2], dofs] = tangents[j] + end + end + + # remove inactive nodes from assembly + for j in S + dofs = [2*(j-1)+1, 2*(j-1)+2] + 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 + + return + +end diff --git a/src/problems_mortar_3d_autodiff.jl b/src/problems_mortar_3d_autodiff.jl deleted file mode 100644 index 052e9cc..0000000 --- a/src/problems_mortar_3d_autodiff.jl +++ /dev/null @@ -1,513 +0,0 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md - -""" Fast inverse of 3x3 matrix. """ -function inv3(P::Matrix) - n, m = size(P) - @assert n == m == 3 - a, b, c, d, e, f, g, h, i = P - A = e*i - f*h - B = -d*i + f*g - C = d*h - e*g - D = -b*i + c*h - E = a*i - c*g - F = -a*h + b*g - G = b*f - c*e - H = -a*f + c*d - I = a*e - b*d - return 1/(a*A + b*B + c*C)*[A B C; D E F; G H I] -end - -""" Project vertex `p` from element surface to auxiliary plane defined -with centerpoint `x0` and normal direction `n0`. -""" -function project_vertex_to_auxiliary_plane(p::Vector, x0::Vector, n0::Vector) - return p - dot(p-x0, n0)*n0 -end - -""" Project vertex `p` from auxiliary plane (x0, n0) back to element surface. - -This requires solving nonlinear system of equations - -f(α,ξ₁,ξ₂) = Nₖξₖ - αn₀ - p = 0 - -""" -function project_vertex_to_surface{E}(p::Vector, x0::Vector, n0::Vector, - element::Element{E}, x::DVTI, time::Real; max_iterations::Int=10, iter_tol::Float64=1.0e-9) - basis(xi) = get_basis(E, xi) - dbasis(xi) = get_dbasis(E, xi) - f(theta) = basis(theta[1:2])*x - theta[3]*n0 - p - L(theta) = inv3([dbasis(theta[1:2])*x -n0]) -# L2(theta) = inv(ForwardDiff.get_value([dbasis(theta[2:3])*x -n0])) - # FIXME: for some reason forwarddiff gives NaN's here. - theta = zeros(3) - dtheta = zeros(3) - for i=1:max_iterations - dtheta = L(theta) * f(theta) - theta -= dtheta - norm(ForwardDiff.get_value(dtheta)) < iter_tol && return theta[1:2], theta[3] - end - - info("failed to project vertex from auxiliary plane back to surface") - info("element type: $E") - info("element connectivity: $(get_connectivity(element))") - info("auxiliary plane: x0 = $(ForwardDiff.get_value(x0)), n0 = $(ForwardDiff.get_value(n0))") - info("element geometry: $(ForwardDiff.get_value(x.data))") - info("vertex to project: $(ForwardDiff.get_value(p))") - info("parameter vector before giving up: $(ForwardDiff.get_value(theta)')") - info("increment in parameter vector before giving up: $(ForwardDiff.get_value(dtheta)')") - info("norm(dtheta) before giving up: $(ForwardDiff.get_value(norm(dtheta)))") - info("f([0.0, 0.0, 0.0]) = $(ForwardDiff.get_value(f([0.0, 0.0, 0.0]))')") - info("L([0.0, 0.0, 0.0]) = $(ForwardDiff.get_value(L([0.0, 0.0, 0.0])))") - - info("iterations:") - theta = zeros(3) - dtheta = zeros(3) - for i=1:max_iterations - info("iter $i, theta = $(ForwardDiff.get_value(theta)')") - info("f = $(ForwardDiff.get_value(f(theta))')") - info("L = $(ForwardDiff.get_value(L(theta)))") -# info("L2 = $(ForwardDiff.get_value(L2(theta)))") - dtheta = L(theta) * f(theta) - info("dtheta = $(ForwardDiff.get_value(dtheta)')") - theta -= dtheta - end - - error("project_point_to_surface: did not converge in $max_iterations iterations!") -end - -""" Test is q inside sm. -http://bbs.dartmouth.edu/~fangq/MATH/download/source/Determining%20if%20a%20point%20lies%20on%20the%20interior%20of%20a%20polygon.htm -""" -function vertex_inside_polygon(q, P; atol=1.0e-6) - N = length(P) - angle = 0.0 - for i=1:N - A = P[i] - q - B = P[mod(i,N)+1] - q - c = norm(A)*norm(B) - isapprox(c, 0.0; atol=atol) && return true - cosa = dot(A,B)/c - isapprox(cosa, 1.0; atol=atol) && return false - isapprox(cosa, -1.0; atol=atol) && return true - try - angle += acos(cosa) - catch - info("Unable to calculate acos($(ForwardDiff.get_value(cosa))) when determining is a vertex inside polygon.") - info("Polygon is: $(ForwardDiff.get_value(P)) and vertex under consideration is $(ForwardDiff.get_value(q))") - info("Polygon corner point in loop: A=$(ForwardDiff.get_value(A)), B=$(ForwardDiff.get_value(B))") - info("c = ||A||*||B|| = $(ForwardDiff.get_value(c))") - rethrow() - end - end - return isapprox(angle, 2*pi; atol=atol) -end - -function calculate_centroid(P) - N = length(P) - P0 = P[1] - areas = [norm(1/2*cross(P[i]-P0, P[mod(i,N)+1]-P0)) for i=2:N] - centroids = [1/3*(P0+P[i]+P[mod(i,N)+1]) for i=2:N] - C = 1/sum(areas)*sum(areas.*centroids) - return C -end - -function get_polygon_clip(xs, xm, n) - # objective: search does line xm1 - xm2 clip xs - nm = length(xm) - ns = length(xs) - P = [] - - # 1. test is master point inside slave, if yes, add to clip - for i=1:nm - vertex_inside_polygon(xm[i], xs) && push!(P, xm[i]) - end - - # 2. test is slave point inside master, if yes, add to clip - for i=1:ns - vertex_inside_polygon(xs[i], xm) && push!(P, xs[i]) - end - - for i=1:nm - # 2. find possible intersection - xm1 = xm[i] - xm2 = xm[mod(i,nm)+1] - #info("intersecting line $xm1 -> $xm2") - for j=1:ns - xs1 = xs[j] - xs2 = xs[mod(j,ns)+1] - #info("clipping polygon edge $xs1 -> $xs2") - tnom = dot(cross(xm1-xs1, xm2-xm1), n) - tdenom = dot(cross(xs2-xs1, xm2-xm1), n) - isapprox(tdenom, 0) && continue - t = tnom/tdenom - (0 <= t <= 1) || continue - q = xs1 + t*(xs2 - xs1) - #info("t=$t, q=$q, q ∈ xm ? $(vertex_inside_polygon(q, xm))") - vertex_inside_polygon(q, xm) && push!(P, q) - end - end - - return P -end - -""" Divide polygon to cells. """ -function get_cells(P, C) - N = length(P) - cells = Vector[] - # shared edge etc. - N < 3 && return cells - # trivial case, polygon already triangle / quadrangle - #N == 3 && return Vector[P] - #N == 4 && return Vector[P] - #V = sum([cross(P[i], P[mod(i,N)+1]) for i=1:N]) - #A = 1/2*abs(dot(n, V)) - #info("A = $A") - cells = Vector[Vector[C, P[i], P[mod(i,N)+1]] for i=1:N] - return cells - - maxa = 0.0 - maxj = 0 - for i=1:N - A = P[i] - C - B = P[mod(i,N)+1] - C - theta = acos(dot(A,B)/(norm(A)*norm(B))) - if theta > maxa - maxa = theta - maxj = i - end - end - info("max angle $(maxa/pi*180) at index $maxj, N=$N") - indices = mod(collect(maxj:maxj+N), N) - info("indices = $indices") -end - -""" Check that polygon P is in CCW order for the direction n. Reorder if not. """ -function check_orientation!(P, n) - C = mean(P) - np = length(P) - s = [dot(n, cross(P[i]-C, P[mod(i+1,np)+1]-C)) for i=1:np] - all(s .< 0) && return - info("polygon not in ccw order, fixing") - # project points to new orthogonal basis Q and sort there - t1 = (P[1]-C)/norm(P[1]-C) - t2 = cross(n, t1) - Q = [n t1 t2] - sort!(P, lt=(A, B) -> begin - A_proj = Q'*(A-C) - B_proj = Q'*(B-C) - a = atan2(A_proj[3], A_proj[2]) - b = atan2(B_proj[3], B_proj[2]) - return a > b - end) -end - -""" Assemble Mortar problem for three-dimensional problems, i.e. for Tri3, Tri6, Quad4, Quad8, Quad9 elements. """ -function assemble!(problem::Problem{Mortar}, time::Real, ::Type{Val{3}}) - - 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) # x = [u; la] - 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) - gap = zeros(u) - C = zeros(la) - all_slave_nodes = Set{Int64}() - slave_surface_area = 0.0 - slave_surface_area_2 = 0.0 - slave_element_areas = [] - - # 1. calculate and average node normals for slave element nodes - normal = zeros(u) - for element in get_elements(problem) - haskey(element, "master elements") || continue - conn = get_connectivity(element) - push!(all_slave_nodes, conn...) - 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) - J = transpose(sum([kron(dN[:,i], x_el[i]') for i=1:length(x_el)])) - n = reshape(cross(J[:,1], J[:,2]), 3, 1) - normal[:, conn] += ip.weight*n*N - end - end - - all_slave_nodes = sort(collect(all_slave_nodes)) - - # normalize to unit normal - for i in all_slave_nodes - normal[:,i] /= norm(normal[:,i]) - end - #normal = ForwardDiff.get_value(normal) - - if props.rotate_normals - for i=1:size(normal, 2) - normal[:,i] = -normal[:,i] - end - end - - # 2. loop slave elements and find contact segments - for slave_element in get_elements(problem) - slave_element_area = 0.0 - haskey(slave_element, "master elements") || continue - info("new slave element") - - slave_element_nodes = get_connectivity(slave_element) - X1 = slave_element("geometry", time) - u1 = Field(Vector[u[:,i] for i in slave_element_nodes]) - if haskey(slave_element, "displacement") - u1 -= slave_element("displacement", time) - end - x1 = X1 + u1 - la1 = Field(Vector[la[:,i] for i in slave_element_nodes]) - #if haskey(slave_element, "reaction force") - # la1 -= slave_element("reaction force", time) - #end - n1 = Field(Vector[normal[:,i] for i in slave_element_nodes]) - nnodes = size(slave_element, 2) - update!(slave_element, "normals", time => ForwardDiff.get_value(n1.data)) - - # 2.1. create auxiliary plane (x0, Q) - xi = get_reference_element_midpoint(slave_element) - N = vec(get_basis(slave_element, xi)) - x0 = N*x1 - n0 = N*n1 - - # 2.2. project slave nodes to auxiliary plane - 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"] - - master_element_nodes = get_connectivity(master_element) - X2 = master_element("geometry", time) - u2 = Field(Vector[u[:,i] for i in master_element_nodes]) - if haskey(master_element, "displacement") - u2 -= master_element("displacement", time) - end - x2 = X2 + u2 - - distance = norm(mean(x2) - mean(x1)) - distance > props.maximum_distance && continue - - # 3.1. project master nodes to auxiliary plane - M = Vector[project_vertex_to_auxiliary_plane(p, x0, n0) for p in x2] - - # 3.2. create polygon clipping on auxiliary plane - 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) - - # 3.3. loop integration cells one at time - for cell in get_cells(P, C0) - x_cell = Field(cell) - - # 3.3.1. create dual basis - De = zeros(nnodes, nnodes) - Me = zeros(nnodes, nnodes) - for ip in get_integration_points(Tri3, Val{5}) - N = vec(get_basis(Tri3, ip.xi)) - x_gauss = N*x_cell - xi_slave, alpha = project_vertex_to_surface(x_gauss, x0, n0, slave_element, x1, time) - N1 = slave_element(xi_slave, time) - dNC = get_dbasis(Tri3, ip.xi) - JC = transpose(sum([kron(dNC[:,j], x_cell[j]') for j=1:length(x_cell)])) - wC = ip.weight*norm(cross(JC[:,1], JC[:,2])) - De += wC*diagm(vec(N1)) - Me += wC*N1'*N1 - end - Ae = De*inv(Me) - - # 3.3.2 loop integration points of cell and calculate fc and gap - for ip in get_integration_points(Tri3, Val{5}) - N = vec(get_basis(Tri3, ip.xi)) - x_gauss = N*x_cell - - # project gauss point back to element surfaces - xi_slave, alpha = project_vertex_to_surface(x_gauss, x0, n0, slave_element, x1, time) - xi_master, alpha = project_vertex_to_surface(x_gauss, x0, n0, master_element, x2, time) - - # evaluate shape functions, calculate contact force and gap - N1 = vec(get_basis(slave_element, xi_slave)) - N2 = vec(get_basis(master_element, xi_master)) - Phi = Ae*N1 - - dNC = get_dbasis(Tri3, ip.xi) - JC = transpose(sum([kron(dNC[:,j], x_cell[j]') for j=1:length(x_cell)])) - wC = ip.weight*norm(cross(JC[:,1], JC[:,2])) - - x_s = N1*x1 - x_m = N2*x2 - u_s = N1*u1 - u_m = N2*u2 - n_s = N1*n1 - la_s = Phi*la1 - la_n = dot(n_s, la_s) - g_s = x_s-x_m - #gn = -dot(n_s, g_s) - fc[:,slave_element_nodes] += wC*la_s*N1' - fc[:,master_element_nodes] -= wC*la_s*N2' - #gap[:,slave_element_nodes] += wC*g_s*N1' - #gap[:,master_element_nodes] += wC*g_s*N2' - gn = props.gap_sign*dot(n_s, g_s) - #gap[1,slave_element_nodes] += wC*gn*Phi' - #gap[1,slave_element_nodes] += wC*gn*Phi' - #C[:,master_element_nodes] -= wC*u_s*N2' - gap[:,slave_element_nodes] = wC*props.gap_sign*g_s*Phi' - #gap[:,master_element_nodes] -= wC*(u_s-u_m)*N2' - - slave_surface_area += wC - slave_element_area += wC - end # done integrating cell - - end # done for all cells in this segment - - end # done all master elements for this slave element - push!(slave_element_areas, slave_element_area) - - end # done all slave elements - - # like in 2d, check contact in nodes based on a complementarity condition - - nzgap = sort(nonzeros(sparse(ForwardDiff.get_value(gap)))) - info("gap: $nzgap") - - info("size of normal = $(size(normal))") - info("size of la = $(size(la))") - info("size of C = $(size(C))") - info("S = $all_slave_nodes") - info("slave surface area: $(ForwardDiff.get_value(slave_surface_area))") - info("slave surface area 2: $(ForwardDiff.get_value(slave_surface_area_2))") - info("slave element areas:") - for (i, a) in enumerate(slave_element_areas) - info("element $i, area = $(ForwardDiff.get_value(a))") - end - for (i, element) in enumerate(get_elements(problem)) - haskey(element, "master elements") || continue - info("element $i geometry: $(element("geometry", time).data)") - end - - for (i, j) in enumerate(all_slave_nodes) - if j in props.always_inactive - info("special node $j always inactive") - C[:,j] = la[:,j] - continue - end - n = normal[:,j] - 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) - Q = [n t1 t2] - la_nt = Q'*la[:,j] - gap_nt = Q'*gap[:,j] - C[1,j] = gap_nt[1] - C[2:3,j] = la_nt[2:3] - #C[:,j] -= gap[:,j] -#= - if lan - gn < 0 - info("set node $j active, normal direction = $(ForwardDiff.get_value(n)), tangent plane = $(ForwardDiff.get_value(t1)) x $(ForwardDiff.get_value(t2))") - C[1,j] = gn - C[2,j] = dot(t1, la[:,j]) - C[3,j] = dot(t2, la[:,j]) - else - C[:,j] = la[:,j] - end -=# - - end - -#= - for (i, j) in enumerate(all_slave_nodes) - n = normal[:,j] - 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) - Q = [n t1 t2] - Ci = ForwardDiff.get_value(Q'*C[:,j]) - gapi = ForwardDiff.get_value(Q'*gap[:,j]) - fci = ForwardDiff.get_value(Q'*fc[:,j]) - lai = ForwardDiff.get_value(Q'*la[:,j]) - ui = ForwardDiff.get_value(Q'*u[:,j]) - info("$i/$j: \nC = $Ci, \nf = $fci, \ngap = $gapi, \nla = $lai, \nu = $ui") - 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) - A, allresults = ForwardDiff.jacobian(calculate_interface, x, - ForwardDiff.AllResults, cache=autodiffcache) - b = ForwardDiff.value(allresults) - -# dump(round(A, 3)) -# dump(round(b, 3)') - A = sparse(A) - b = sparse(b) - - SparseMatrix.droptol!(A, 1.0e-12) - SparseMatrix.droptol!(b, 1.0e-12) - - 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] - - function joo(x) - nz1 = sort(unique(rowvals(x))) - nz2 = sort(unique(rowvals(x'))) - info("nz1 = $nz1, nz2 = $nz2") - dump(round(full(x[nz1,nz2]), 3)) - end - println("K") - joo(K) - println("C1") - joo(C1) - println("C2") - joo(C2) - println("D") - joo(D) - println("f") - joo(f) - println("g") - joo(g) -#= - slaves = [101,108,111,112,113,120,123,124,125,126,129,130,149,150,151,152] - for j in slaves - dofs = [3*(j-1)+1, 3*(j-1)+2, 3*(j-1)+3] - info("slave node $j, dofs $dofs") - info("Stiffness: $(K[dofs,:])") - info("force fc: $(C1[dofs,:])") - info("constraint: $(C2[dofs,:])") - info("lambdas: $(D[dofs,:])") - info("f = $(f[dofs]), g = $(g[dofs])") - end -=# - - empty!(problem.assembly) - add!(problem.assembly.K, K) - add!(problem.assembly.C1, C1) - add!(problem.assembly.C2, C2) - add!(problem.assembly.D, D) - add!(problem.assembly.f, f) - add!(problem.assembly.g, g) - - return problem.assembly -end diff --git a/src/problems_mortar_forwarddiff_pe.jl b/src/problems_mortar_forwarddiff_pe.jl deleted file mode 100644 index f0117a8..0000000 --- a/src/problems_mortar_forwarddiff_pe.jl +++ /dev/null @@ -1,446 +0,0 @@ -using JuliaFEM.Core: MortarElements2D, DVTI, Assembly - -import JuliaFEM.Core: project_from_master_to_slave, project_from_slave_to_master, assemble!, -get_unknown_field_dimension, get_parent_field_name, get_gdofs, find_elements, get_nodes, Field, -get_integration_points, get_basis, get_dbasis, add! - -""" Find segment from slave element corresponding to master element nodes. -x1_, n1_ -slave element geometry and normal direction - -x2_ master element nodes to project onto slave -""" -function project_from_master_to_slave{E<:MortarElements2D}( - slave_element::Element{E}, x1_::DVTI, n1_::DVTI, x2::Vector) - - function x1(xi1) - N = get_basis(E, xi1) - return vec(N)*x1_ - end - - function dx1(xi1) - dN = get_dbasis(E, xi1) - return vec(dN)*x1_ - end - - function n1(xi1) - N = get_basis(E, xi1) - return vec(N)*n1_ - end - - function dn1(xi1) - dN = get_dbasis(E, xi1) - return vec(dN)*n1_ - end - - cross2(a, b) = cross([a; 0], [b; 0])[3] - R(xi1) = cross2(x1(xi1)-x2, n1(xi1)) - dR(xi1) = cross2(dx1(xi1), n1(xi1)) + cross2(x1(xi1)-x2, dn1(xi1)) - xi1 = 0.0 - for i=1:5 - dxi1 = -R(xi1)/dR(xi1) - xi1 += dxi1 - if norm(dxi1) < 1.0e-10 - return xi1 - end - end - - error("find projection from master to slave: did not converge") - -end - -function project_from_slave_to_master{E<:MortarElements2D}( - master_element::Element{E}, x1::Vector, n1::Vector, x2_::DVTI) - - function x2(xi2) - N = get_basis(E, xi2) - return vec(N)*x2_ - end - - function dx2(xi2) - dN = get_dbasis(E, xi2) - return vec(dN)*x2_ - end - - cross2(a, b) = cross([a; 0], [b; 0])[3] - R(xi2) = cross2(x2(xi2)-x1, n1) - dR(xi2) = cross2(dx2(xi2), n1) - xi2 = 0.0 - dxi2 = 0.0 - for i=1:5 - dxi2 = -R(xi2) / dR(xi2) - xi2 += dxi2 - if norm(dxi2) < 1.0e-10 - return xi2 - end - 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}}) - 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(u::Matrix, la::Matrix) - - X1 = slave_element("geometry", time) - slave_element_nodes = get_connectivity(slave_element) - 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 - - 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] - # 1. update nodal normals for this element - 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]) - - fc = SparseMatrixCOO{Real}([], [], []) # interface virtual work - C = SparseMatrixCOO{Real}([], [], []) # constraints - B = SparseMatrixCOO{Real}([], [], []) - - #info("u1.data = ", ForwardDiff.get_value(u1.data)) - info("normal calculations done. looping master elements.") - for master_element in slave_element["master elements"] - X2 = master_element("geometry", time) - master_element_nodes = get_connectivity(master_element) - u2 = Field(Vector[u[:,i] for i in master_element_nodes]) - x2 = X2 + u2 - info("master element ready.") - - # calculate segmentation: we care only about endpoints - # note: these are quadratic/cubic functions, analytical solution possible - info("calculating segmentation.") - 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 - info("xi1 = $xi1") - - info("create bi-orthogonal basis") - 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) - info("bi-orthogonal basis done. integrating fc.") - - slave_dofs = get_gdofs(slave_element, field_dim) - master_dofs = get_gdofs(master_element, field_dim) - - info("integrate fc") - D = zeros(nnodes, nnodes) - M = zeros(nnodes, nnodes) - gn = zeros(nnodes) - lan = zeros(nnodes) - lat = zeros(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) - # project gauss point to master element to evaluate shape function there - x_s = vec(N1)*x1 # coordinate in gauss point - n_s = vec(N1)*n1 # normal direction in gauss point - xi_m = project_from_slave_to_master(master_element, x_s, n_s, x2) - N2 = get_basis(master_element, xi_m) - x_m = vec(N2)*x2 - Phi = vec(Ae*N1') - la_s = Phi*la1 # traction force in gauss point - u_s = vec(N1)*u1 - u_m = vec(N2)*u2 - #info("la_s = $(ForwardDiff.get_value(la_s))") - SM = [slave_dofs; master_dofs] - #N1N2 = [N1 -N2] - #info("all_dofs = $(SM)") - #info("shape functions = $(ForwardDiff.get_value(N1N2))") - #for i=1:field_dim - # #info("add to slave dofs $(slave_dofs[i:field_dim:end])") - # #info("add to master dofs $(master_dofs[i:field_dim:end])") - # add!(fc, slave_dofs[i:field_dim:end], [1, 1], -w*la_s[i]*N1) - # add!(fc, master_dofs[i:field_dim:end], [1, 1], +w*la_s[i]*N2) - #add!(fc, slave_dofs[i:field_dim:end], [1, 1], w*la_s'*u_s[i]) - #add!(fc, master_dofs[i:field_dim:end], [1, 1], -w*la_s'*u_m[i]) - # add!(fc, master_dofs[i:field_dim:end], [1, 1], -w*la_s[i]*u_m) - #end - #add!(fc, [slave_dofs; master_dofs], [1, 1, 1, 1, 1, 1, 1, 1], w*la_s*[u_s' -u_m']) - D += w*kron(Ae*N1', N1) - M += w*kron(Ae*N1', N2) - gn += -w*dot(n_s, x_s-x_m)*Phi - lan += w*dot(n_s, la_s)*Phi - t_s = Q'*n_s - lat += w*dot(t_s, la_s)*Phi - end - - #D2 = zeros(2*nnodes, 2*nnodes) - #M2 = zeros(2*nnodes, 2*nnodes) - #for i=1:field_dim - # D2[i:field_dim:end, i:field_dim:end] += D - # M2[i:field_dim:end, i:field_dim:end] += M - #end - #info("size of D2 = $(size(D2))") - #fco = [D2 -M2]*vec(la1) - #fco = [D -M]*la[:,slave_element_nodes] - #info("fco = $(ForwardDiff.get_value(fco))") - #add!(fc, [slave_dofs; master_dofs], [1, 1, 1, 1], fco) - - 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 - info("gn = $gn") - #Cj = lan - max(0, lan - gn) + lat - add!(C, slave_dofs[1:field_dim:end], [1, 1], gn') - - end # master elements done - - ndofs = prod(size(la)) - 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[1]], dofs, reshape(n1[i], 1, 2)) - add!(T, [dofs[2]], dofs, reshape(Q'*n1[i], 1, 2)) - end - N = sparse(N, ndofs, ndofs) - T = sparse(T, ndofs, ndofs) - B = sparse(B, ndofs, ndofs) - fc = B'*vec(la) - #println(sparse(fc)) - #fc = sparse(fc, ndofs, 1) - #println(fc) - #dump(full(fc)) - #C = sparse(C, ndofs, 1) - C = N*B*vec(u) + T*vec(la) - return fc, C - - end - - - function calculate_interface_PE(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) - #fixed_la = ForwardDiff.get_value(la) - #fixed_u = ForwardDiff.get_value(u) - #u = ForwardDiff.get_value(u) - - X1 = slave_element("geometry", time) - slave_element_nodes = get_connectivity(slave_element) - u1 = Field(Vector[u[:,i] for i in slave_element_nodes]) - la1 = Field(Vector[la[:,i] for i in slave_element_nodes]) - #fixed_la1 = Field(Vector[fixed_la[:,i] for i in slave_element_nodes]) - x1 = X1 + u1 - - # 1. update nodal normals for this 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] - 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]) - - Wco = 0.0 - Wla = 0.0 - - for master_element in slave_element["master elements"] - X2 = master_element("geometry", time) - master_element_nodes = get_connectivity(master_element) - 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 - - 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) - - 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 = vec(get_basis(slave_element, xi_s)) - - # project gauss point to master element to evaluate shape function there - x_s = N1*x1 # coordinate in gauss point - n_s = N1*n1 # normal direction in gauss point - t_s = Q'*n_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 - - gn = -dot(n_s, x_s - x_m) - gt = dot(t_s, x_s - x_m) - lan = dot(n_s, Phi*la1) - lat = dot(t_s, Phi*la1) - u_s = N1*u1 - u_m = N2*u2 - gu = dot(n_s, u_s - u_m) - - Wco += w*dot(Phi*la1, N1*u1 - N2*u2) - #Wla += w*(lan*gn + lat*gt) - #gn = min(0, gn) - Wla += 1/2*w*1e6*gn*gn - #info("gn = $(ForwardDiff.get_value(gn))") - #Wla += w*dot(dot(n_s, Phi*la1), dot(n_s, N1*u1 - N2*u2)) - end - - end - return Wco, Wla - end - - - function calculate_contact_rhs(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, C = calculate_interface(u, la) - info("interface vector calculated.") - return vec(full([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 - - function add_fco!() - get_PI(x::Vector) = calculate_interface_PE(x)[1] - A, allresults = ForwardDiff.hessian(get_PI, x, ForwardDiff.AllResults) - b = -ForwardDiff.gradient(allresults) - info("PE = $(ForwardDiff.value(allresults))") - A = sparse(A) - b = sparse(b) - SparseMatrix.droptol!(A, 1.0e-12) - SparseMatrix.droptol!(b, 1.0e-12) - 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) - end - #add_fco!() - - function add_wla!() - get_PI(x::Vector) = calculate_interface_PE(x)[2] - A, allresults = ForwardDiff.hessian(get_PI, x, ForwardDiff.AllResults) - b = -ForwardDiff.gradient(allresults) - info("PE = $(ForwardDiff.value(allresults))") - A = sparse(A) - b = sparse(b) - SparseMatrix.droptol!(A, 1.0e-12) - SparseMatrix.droptol!(b, 1.0e-12) - #info("A") - #println(full(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) - end - add_wla!() - - return - -end - -mesh, body1, body2, bc_top, bc_bottom, contact = divided_block_problem() -bc_top.properties.formulation = :incremental -bc_bottom.properties.formulation = :incremental -contact.properties.formulation = :forwarddiff -contact.assembly.u = zeros(16) -contact.assembly.la = zeros(16) -assemble!(contact.assembly, contact, contact.elements[1], 0.0, Val{:forwarddiff}) -