diff --git a/src/mortar.jl b/src/mortar.jl index e1c8246..9e3999b 100644 --- a/src/mortar.jl +++ b/src/mortar.jl @@ -24,18 +24,20 @@ type Mortar <: BoundaryProblem inequality_constraints :: Bool # Launch PDASS to solve inequality constraints normal_condition :: Symbol # Tie or Contact tangential_condition :: Symbol # Stick or Slip - minimum_distance :: Float64 # don't check for a contact if elements are far enough + maximum_distance :: Float64 # don't check for a contact if elements are far enough store_debug_info :: Bool # for making debugging easier + always_inactive :: Vector{Int64} always_in_contact :: Vector{Int64} # nodes in this list always in contact always_in_stick :: Vector{Int64} # nodes in this list always in stick always_in_slip :: Vector{Int64} # nodes in this list always in slip contact :: Bool friction :: Bool gap_sign :: Int # gap sign convention + rotate_normals :: Bool end function Mortar() - Mortar(:total, true, false, :Tie, :Stick, Inf, false, [], [], [], false, false, -1) + Mortar(:total, true, false, :Tie, :Stick, Inf, false, [], [], [], [], false, false, -1, false) end function get_unknown_field_name(::Type{Mortar}) diff --git a/src/mortar_2d_autodiff.jl b/src/mortar_2d_autodiff.jl index f638f8f..d8ea24b 100644 --- a/src/mortar_2d_autodiff.jl +++ b/src/mortar_2d_autodiff.jl @@ -73,10 +73,8 @@ function project_from_slave_to_master{E<:MortarElements2D}( 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 +function assemble!(problem::Problem{Mortar}, time::Real) + props = problem.properties field_dim = get_unknown_field_dimension(problem) field_name = get_parent_field_name(problem) @@ -88,21 +86,17 @@ function assemble!{E<:MortarElements2D}(assembly::Assembly, 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) + S = Set{Int64}() - 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 + # 1. update nodal normals for slave elements Q = [0.0 -1.0; 1.0 0.0] normals = zeros(u) - for element in adjacent_elements + for element in get_elements(problem) + haskey(element, "master elements") || continue conn = get_connectivity(element) + push!(S, conn...) gdofs = get_gdofs(element, field_dim) X_el = element("geometry", time) u_el = Field(Vector[u[:, i] for i in conn]) @@ -114,94 +108,129 @@ function assemble!{E<:MortarElements2D}(assembly::Assembly, 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]) - - nnodes = size(slave_element, 2) - lan_tot = zeros(nnodes) # normal pressure - gap_tot = zeros(nnodes) # weighted normal gap - - 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 - - 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 + for i in 1:size(normals,2) + normals[:,i] /= norm(normals[:,i]) + end + if props.rotate_normals + for i=1:size(normals,2) + normals[:,i] = -normals[:,i] end - Ae = De*inv(Me) + end - slave_dofs = get_gdofs(slave_element, field_dim) - master_dofs = get_gdofs(master_element, field_dim) + # 2. loop all slave elements + for slave_element in get_elements(problem) + haskey(slave_element, "master elements") || continue - 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 + slave_element_nodes = get_connectivity(slave_element) + X1 = slave_element("geometry", time) + u1 = Field(Vector[u[:,i] for i in slave_element_nodes]) + x1 = X1 + u1 + la1 = Field(Vector[la[:,i] for i in slave_element_nodes]) + n1 = Field(Vector[normals[:,i] for i in slave_element_nodes]) - # 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 + nnodes = size(slave_element, 2) - u_s = N1*u1 - u_m = N2*u2 + # 3. loop all master elements + for master_element in slave_element["master elements"] - la_s = Phi*la1 # traction force in gauss point - la_nt = R_s*la_s - gn = -dot(n_s, x_s - x_m) # normal gap + 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 - fc[:,slave_element_nodes] += w*la_s*N1' - fc[:,master_element_nodes] -= w*la_s*N2' - #C[1,slave_element_nodes] += w*gn*Phi' + x1_midpoint = 1/2*(x1[1]+x1[2]) + x2_midpoint = 1/2*(x2[1]+x2[2]) + distance = ForwardDiff.get_value(norm(x2_midpoint - x1_midpoint)) + distance > props.maximum_distance && continue - lan_tot += w*la_nt[1]*Phi - gap_tot += w*gn*Phi + # calculate segmentation: we care only about endpoints + # note: these are quadratic/cubic functions, analytical solution possible + xi1a = -Inf + xi1b = -Inf + try + xi1a = project_from_master_to_slave(slave_element, x1, n1, x2[1]) + xi1b = project_from_master_to_slave(slave_element, x1, n1, x2[end]) + catch + info("failed to create projection!!!!") + # TODO + continue + 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 + 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) + + # 4. loop integration points of segment + 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 + 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 + + la_s = Phi*la1 # traction force in gauss point + gn = props.gap_sign*dot(n_s, x_s - x_m) # normal gap + + fc[:,slave_element_nodes] += w*la_s*N1' + fc[:,master_element_nodes] -= w*la_s*N2' + gap[1,slave_element_nodes] += w*gn*Phi' + #gap[1,slave_element_nodes] += w*gn*N1' + + end + + end # master elements done + + end # slave elements done + + # at this point we have calculated contact force fc and gap for all slave elements. + # next task is to find out are they in contact or not and remove inactive nodes + + nzgap = sort(nonzeros(sparse(ForwardDiff.get_value(gap)))) + info("gap: $nzgap") + + for (i, j) in enumerate(sort(collect(S))) + if j in props.always_inactive + info("special node $j always inactive") + C[:,j] = la[:,j] + continue end - - end # master elements done - - for (i, j) in enumerate(slave_element_nodes) - n = n1[i] + n = normals[:,j] t = Q'*n - R = [n t] - la_nt = R'*la[:,j] -# info("node $j, n=$(ForwardDiff.get_value(n)) lan = $(ForwardDiff.get_value(la_nt[1])) gap = $(ForwardDiff.get_value(gap_tot[i]))") + lan = dot(n, la[:,j]) + lat = dot(t, la[:,j]) - if la_nt[1] - gap_tot[i] > 0 - info("set node $j active") - C[1,j] += gap_tot[i] - C[2,j] += la_nt[2] + if lan - gap[1, j] > 0 + info("set node $j active, normal direction = $(ForwardDiff.get_value(n)), tangent plane = $(ForwardDiff.get_value(t))") + C[1,j] += gap[1, j] + C[2,j] += lat else - info("set node $j inactive") + #info("set node $j inactive") C[:,j] = la[:,j] end end @@ -213,28 +242,30 @@ function assemble!{E<:MortarElements2D}(assembly::Assembly, # 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) + A, allresults = ForwardDiff.jacobian(calculate_interface, x, + ForwardDiff.AllResults, cache=autodiffcache) 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 + 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/solver_utils.jl b/src/solver_utils.jl index 61ad64f..4987fa8 100644 --- a/src/solver_utils.jl +++ b/src/solver_utils.jl @@ -70,6 +70,7 @@ function handle_overconstraint_error!(problem, nodes, all_dofs, C1_, C1, C2_, C2 # INFO: fixed: new setting is # INFO: dof 1109: 0.0*u₁₅₃ - 0.0*u₁₅₄ - 0.0*u₁₅₅ + 0.15*u₁₅₆ + 0.0*u₁₁₀₉ - 0.15*u₁₁₁₀ = -0.0 =# +#= if 555 in nodes info("overconstraint DIRTY HACK") # It is possible to selectively remove mortar constraints and the associated @@ -103,6 +104,7 @@ function handle_overconstraint_error!(problem, nodes, all_dofs, C1_, C1, C2_, C2 #C1_[1110,:] = 0 return end +=# """ Return all other dofs which connects to overconstrained dofs. """ function get_related_dofs(dofs) @@ -243,7 +245,7 @@ function handle_overconstraint_error!(problem, nodes, all_dofs, C1_, C1, C2_, C2 return dofs_, true end - actions = [action1, action2] + actions = [action1, action2, action3, action4] function show_lambda_coefficients(dofs, C1) for dof in dofs