From 0d0d200bfc36be8147bb2d9055301d883ae36014 Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Tue, 16 Feb 2016 17:40:11 +0200 Subject: [PATCH] Incremental formulation works now, also finite sliding without proper linearization (very slow convergence) --- src/dirichlet.jl | 35 +++++--- src/elements.jl | 6 +- src/mortar.jl | 11 +-- src/mortar_2d.jl | 215 ++++++++++++++++++++++++++++++++++++++++++++++- src/problems.jl | 11 +-- src/solvers.jl | 12 ++- 6 files changed, 256 insertions(+), 34 deletions(-) diff --git a/src/dirichlet.jl b/src/dirichlet.jl index 06c84b1..333726c 100644 --- a/src/dirichlet.jl +++ b/src/dirichlet.jl @@ -1,19 +1,27 @@ # This file is a part of JuliaFEM. # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md +""" Here formulation is :total or :incremental meaning that we either give +constraint for total quantity u or it's increment Δu. For elasticity we are +using incremental formulation. +""" type Dirichlet <: BoundaryProblem formulation :: Symbol dual_basis :: Bool end function Dirichlet() - Dirichlet(:Equality, true) + Dirichlet(:total, true) end function get_unknown_field_name(::Type{Dirichlet}) return "reaction force" end +function get_formulation_type(problem::Problem{Dirichlet}) + return problem.properties.formulation +end + function assemble!(assembly::Assembly, problem::Problem{Dirichlet}, element::Element, time::Real) @assert problem.properties.dual_basis @@ -28,7 +36,7 @@ function assemble!(assembly::Assembly, problem::Problem{Dirichlet}, element::Ele # left hand side for i=1:field_dim ldofs = gdofs[i:field_dim:end] - if haskey(element, field_name*" $i") || haskey(element, field_name) + if haskey(element, field_name*" $i") add!(assembly.C1, ldofs, ldofs, De) add!(assembly.C2, ldofs, ldofs, De) end @@ -47,21 +55,20 @@ function assemble!(assembly::Assembly, problem::Problem{Dirichlet}, element::Ele N = element(ip, time) Phi = (Ae*N')' - if haskey(element, field_name) - for i=1:field_dim - g = element(field_name, ip, time) - ldofs = gdofs[i:field_dim:end] + g_prev = element(field_name, ip, time) + #info("g_prev = $g_prev") + for i=1:field_dim + ldofs = gdofs[i:field_dim:end] + if haskey(element, field_name*" $i") + g = element(field_name*" $i", ip, time) + if get_formulation_type(problem) == :incremental + g = g - g_prev[i] + end + #info("g_new = $g") add!(assembly.g, ldofs, w*g*Phi') end - else - for i=1:field_dim - ldofs = gdofs[i:field_dim:end] - if haskey(element, field_name*" $i") - g = element(field_name*" $i", ip, time) - add!(assembly.g, ldofs, w*g*Phi') - end - end end + end end diff --git a/src/elements.jl b/src/elements.jl index bf02add..608cfc8 100644 --- a/src/elements.jl +++ b/src/elements.jl @@ -268,7 +268,7 @@ function get_nodes(elements::Vector) return nodes end -""" Calculate normal-tangential coordinates for a set of elements. +""" Calculate normal-tangential coordinates for a set of elements. Notes ----- @@ -314,7 +314,7 @@ function calculate_normal_tangential_coordinates!(elements::Vector, time::Real, node_ids = get_connectivity(element) Q = Matrix{Float64}[ [n[:,i] t[:,i]] for i in node_ids] element["normal-tangential coordinates"] = (time => Q) - element["normals"] = Vector{Float64}[n[:,i] for i in node_ids] + element["normals"] = (time => Vector{Float64}[n[:,i] for i in node_ids]) end end @@ -350,7 +350,7 @@ function calculate_normal_tangential_coordinates!(elements::Vector, time::Real, node_ids = get_connectivity(element) Q = Matrix{Float64}[ [n[:,i] t1[:,i] t2[:,i]] for i in node_ids] element["normal-tangential coordinates"] = (time => Q) - element["normals"] = Vector{Float64}[n[:,i] for i in node_ids] + element["normals"] = (time => Vector{Float64}[n[:,i] for i in node_ids]) end end diff --git a/src/mortar.jl b/src/mortar.jl index c8e82f2..c2aee87 100644 --- a/src/mortar.jl +++ b/src/mortar.jl @@ -19,7 +19,8 @@ b) Remove inactive inequality constraints in assembly level. This is done in """ type Mortar <: BoundaryProblem - formulation :: Symbol # Dual or Standard + formulation :: Symbol # :total or :incremental + dual_basis :: Bool inequality_constraints :: Bool # Launch PDASS to solve inequality constraints normal_condition :: Symbol # Tie or Contact tangential_condition :: Symbol # Stick or Slip @@ -33,7 +34,7 @@ type Mortar <: BoundaryProblem end function Mortar() - Mortar(:Dual, false, :Tie, :Stick, Inf, false, [], [], [], false, false) + Mortar(:total, true, false, :Tie, :Stick, Inf, false, [], [], [], false, false) end function get_unknown_field_name(::Type{Mortar}) @@ -41,7 +42,7 @@ function get_unknown_field_name(::Type{Mortar}) end function get_formulation_type(problem::Problem{Mortar}) - return :incremental + return problem.properties.formulation end macro debug(msg) @@ -686,12 +687,12 @@ function assemble!{E<:MortarElements3D}(assembly::Assembly, problem::Problem{Mor # extend matrices according to the problem dimension (3) @assert length(slave_dofs) == length(master_dofs) Me = wC*Ae*N1'*N2 - for k=1:field_dim + for k=1:field_dim C1M3[k:field_dim:end,k:field_dim:end] += Me end end end # integration of mortar matrices done. - + # constraints in normal-tangential direction and initial weighted gap X1 = vec(slave_element("geometry", time)) X2 = vec(master_element("geometry", time)) diff --git a/src/mortar_2d.jl b/src/mortar_2d.jl index 7ae29ae..64845eb 100644 --- a/src/mortar_2d.jl +++ b/src/mortar_2d.jl @@ -223,10 +223,16 @@ end # quadratic not tested yet typealias MortarElements2D Union{Seg2} +function assemble!{E<:MortarElements2D}(assembly::Assembly, problem::Problem{Mortar}, + slave_element::Element{E}, time::Real) + # for finite deformation we need to use incremental formulation + assemble!(assembly, problem, slave_element, time, Val{problem.properties.formulation}) +end + """ Assemble 2d mortar contribution. Mortar matrices are assembled at initial configuration X, so this works for tie contact and small sliding contact. """ function assemble!{E<:MortarElements2D}(assembly::Assembly, problem::Problem{Mortar}, - slave_element::Element{E}, time::Real) + slave_element::Element{E}, time::Real, ::Type{Val{:total}}) # slave element must have a set of master elements haskey(slave_element, "master elements") || return @@ -292,7 +298,7 @@ function assemble!{E<:MortarElements2D}(assembly::Assembly, problem::Problem{Mor Ae = zeros(nnodes, nnodes) De = zeros(nnodes, nnodes) Me = zeros(nnodes, nnodes) - if problem.properties.formulation == :Dual # Construct dual basis + if problem.properties.dual_basis # Construct dual basis for ip in get_integration_points(slave_element, Val{5}) J = get_jacobian(slave_element, ip, time) w = ip.weight*norm(J)*l @@ -422,6 +428,210 @@ function assemble!{E<:MortarElements2D}(assembly::Assembly, problem::Problem{Mor end +function assemble!{E<:MortarElements2D}(assembly::Assembly, problem::Problem{Mortar}, + slave_element::Element{E}, time::Real, ::Type{Val{:incremental}}) + + # slave element must have a set of master elements + haskey(slave_element, "master elements") || return + props = problem.properties + + # get dimension and name of PARENT field + field_dim = problem.dimension + field_name = problem.parent_field_name + + slave_dofs = get_gdofs(slave_element, field_dim) + nnodes = size(slave_element, 2) + + # slave side quantities: rotation matrix, geometry, displacement, reaction force + Q = slave_element("normal-tangential coordinates", time) + Z = zeros(nnodes, nnodes) + if nnodes == 2 + Q2 = [Q[1] Z; Z Q[2]] + elseif nnodes == 3 + Q2 = [Q[1] Z Z; Z Q[2] Z; Z Z Q[3]] + end + X1 = vec(slave_element("geometry", time)) + u1 = zeros(2*nnodes) + if haskey(slave_element, "displacement") + u1 = vec(slave_element("displacement", time)) + end + x1 = X1 + u1 + la = zeros(2*nnodes) + if haskey(slave_element, "reaction force") + la = vec(slave_element("reaction force", time)) + end + la = Q2'*la + + G = zeros(2*nnodes) + g = zeros(2*nnodes) + local_assembly = Assembly() + + has_contribution = false + + for master_element in slave_element["master elements"] + + X2 = vec(master_element("geometry", time)) + u2 = zeros(2*nnodes) + if haskey(master_element, "displacement") + u2 = vec(master_element("displacement", time)) + end + x2 = X2 + u2 + + # if distance between elements is "far enough" cannot expect contact + if props.contact && (props.minimum_distance < Inf) + slave_midpoint = Float64[mean(x1[1:field_dim:2]), mean(x1[2:field_dim:2])] + master_midpoint = Float64[mean(x2[1:field_dim:2]), mean(x2[2:field_dim:2])] + if norm(slave_midpoint - master_midpoint) > props.minimum_distance + continue + end + end + + master_dofs = get_gdofs(master_element, field_dim) + xi1a = project_from_master_to_slave(slave_element, master_element, [-1.0], time, Val{:deformed}) + xi1b = project_from_master_to_slave(slave_element, master_element, [ 1.0], time, Val{:deformed}) + xi1 = clamp([xi1a xi1b], -1.0, 1.0) + l = 1/2*abs(xi1[2]-xi1[1]) + isapprox(l, 0.0) && continue # no contribution + + # Calculate slave side projection matrix D + Ae = zeros(nnodes, nnodes) + De = zeros(nnodes, nnodes) + Me = zeros(nnodes, nnodes) + if problem.properties.dual_basis # Construct dual basis + for ip in get_integration_points(slave_element, Val{5}) + J = get_jacobian(slave_element, ip, time, Val{:deformed}) + w = ip.weight*norm(J)*l + xi = 1/2*(1-ip.xi)*xi1[1] + 1/2*(1+ip.xi)*xi1[2] + N = slave_element(xi, time) + De += w*diagm(vec(N)) + Me += w*N'*N + end + Ae = De*inv(Me) + else # Standard Lagrange basis + for ip in get_integration_points(slave_element, Val{5}) + J = get_jacobian(slave_element, ip, time, Val{:deformed}) + w = ip.weight*norm(J)*l + xi = 1/2*(1-ip.xi)*xi1[1] + 1/2*(1+ip.xi)*xi1[2] + N = slave_element(xi, time) + De += w*N'*N + end + Ae = eye(nnodes) + end + + C1S2 = zeros(2*nnodes, 2*nnodes) + C1M2 = zeros(2*nnodes, 2*nnodes) + + # Slave side already done; it's De + for i=1:field_dim + C1S2[i:field_dim:end,i:field_dim:end] += De + end + + # Calculate master side projection matrix M + for ip in get_integration_points(slave_element, Val{5}) + J = get_jacobian(slave_element, ip, time, Val{:deformed}) + w = ip.weight*norm(J)*l + # integration point on slave side segment + xi_slave = 1/2*(1-ip.xi)*xi1[1] + 1/2*(1+ip.xi)*xi1[2] + # projected integration point to master side element + xi_master = project_from_slave_to_master(slave_element, master_element, + xi_slave, time, Val{:deformed}) + N1 = slave_element(xi_slave, time) + N2 = master_element(xi_master, time) + M = w*kron(Ae*N1', N2) + for i=1:field_dim + C1M2[i:field_dim:end,i:field_dim:end] += M + end + end + + # Calculate normal-tangential constraints and weighted gap + C2S2 = Q2'*C1S2 + C2M2 = Q2'*C1M2 + G += -(C2S2*X1 - C2M2*X2) + g += -(C2S2*x1 - C2M2*x2) + + # Add contributions + add!(local_assembly.C1, slave_dofs, slave_dofs, C1S2) + add!(local_assembly.C1, slave_dofs, master_dofs, -C1M2) + add!(local_assembly.C2, slave_dofs, slave_dofs, C2S2) + add!(local_assembly.C2, slave_dofs, master_dofs, -C2M2) + has_contribution = true + + end # all master elements are done + + if !has_contribution + return + end + + add!(local_assembly.g, slave_dofs, g) + + # if only equality constraints, i.e., mesh tying problem, we're done for this element. + if !props.contact + append!(assembly, local_assembly) + return + end + + lan = la[1:field_dim:end] + lat = la[2:field_dim:end] + gn = g[1:field_dim:end] + gt = g[2:field_dim:end] + + # normal condition + cn = 1.0 # complemementarity parameter + Cn = lan - max(0, lan - cn*gn) + inactive_nodes = find(lan - cn*gn .<= 0) + active_nodes = find(lan - cn*gn .> 0) + + # if all nodes inactive, nothing to contribute. + if length(active_nodes) == 0 + return + end + + # manipulate local assembly (remove rows from it based on active set) + # before adding it to global assembly + C1 = sparse(local_assembly.C1) + C2 = sparse(local_assembly.C2) + D = spzeros(size(C2)...) + g = sparse(local_assembly.g) + + node_ids = get_connectivity(slave_element) + + # normal constraint: remove inactive nodes + for j in node_ids[inactive_nodes] + if length(props.always_in_contact) != 0 + j in props.always_in_contact && continue + end + gdofs = [2*(j-1)+1, 2*(j-1)+2] + # λⱼ = 0 ∀ j ∈ S + C1[gdofs,:] = 0 + C2[gdofs,:] = 0 + D[gdofs,:] = 0 + g[gdofs,:] = 0 + end + + for (i, j) in enumerate(node_ids[active_nodes]) + gdofs = [2*(j-1)+1, 2*(j-1)+2] + #D[gdofs[2],gdofs] = C2[gdofs[2],gdofs] + D[gdofs[2],gdofs] = Q[i][:,2] + C2[gdofs[2],:] = 0 + g[gdofs[2],:] = 0 + end + + local_assembly.C1 = C1 + local_assembly.C2 = C2 + local_assembly.D = D + local_assembly.g = g + append!(assembly, local_assembly) + + if props.store_debug_info + slave_element["g"] = g + slave_element["c"] = c + slave_element["C1"] = C1 + slave_element["C2"] = C2 + slave_element["D"] = D + slave_element["active nodes"] = active_nodes + end + +end function calculate_gap_vector{E<:MortarElements2D}( problem::Problem{Mortar}, slave_element::Element{E}, @@ -493,4 +703,3 @@ function calculate_gap_vector{E<:MortarElements2D}( return gap end - diff --git a/src/problems.jl b/src/problems.jl index 3caa508..2cd5371 100644 --- a/src/problems.jl +++ b/src/problems.jl @@ -131,7 +131,7 @@ function initialize!(problem::Problem, time::Real) # if this is boundary problem and not dirichlet problem, initialize field # for primary variable too is_boundary_problem(problem) || return - is_dirichlet_problem(problem) && return + #is_dirichlet_problem(problem) && return field_name = get_parent_field_name(problem) for element in get_elements(problem) gdofs = get_gdofs(element, problem) @@ -167,10 +167,12 @@ function update_assembly!(problem, u, la) assembly.u_prev = copy(assembly.u) assembly.la_prev = copy(assembly.la) if get_formulation_type(problem) == :incremental - info("incremental formulation, adding increment to solution vector") + info("$(problem.name): incremental formulation, adding increment to solution vector") + #info("solution vector:") + #dump(round(u, 3)') assembly.u += u else - info("total formulation, replacing solution vector with new values") + info("$(problem.name): total formulation, replacing solution vector with new values") assembly.u = u end assembly.la = la @@ -212,7 +214,7 @@ function update_elements!{P<:BoundaryProblem}(problem::Problem{P}, u, la) last(element[field_name]).data = local_sol end # if boundary problem is not dirichlet, update also data of main problem - is_dirichlet_problem(problem) && return + # is_dirichlet_problem(problem) && return field_name = get_parent_field_name(problem) solution = reshape(u, field_dim, nnodes) for element in get_elements(problem) @@ -285,4 +287,3 @@ function find_nodes_by_dofs(dim, dofs) end return nodes end - diff --git a/src/solvers.jl b/src/solvers.jl index ee24585..22aa4af 100644 --- a/src/solvers.jl +++ b/src/solvers.jl @@ -146,6 +146,7 @@ type Solver is_linear_system :: Bool # setting this to true makes assumption of one step convergence nonlinear_system_max_iterations :: Int64 nonlinear_system_convergence_tolerance :: Float64 + nonlinear_system_error_if_no_convergence :: Bool linear_system_solver :: Symbol end @@ -159,6 +160,7 @@ function Solver(name::ASCIIString="default solver", time::Real=0.0) false, # is_linear_system 10, # max nonlinear iterations 5.0e-5, # nonlinear iteration convergence tolerance + true, # throw error if no convergence :DirectLinearSolver # linear system solution method ) end @@ -344,7 +346,7 @@ function has_converged(solver::Solver; check_convergence_for_boundary_problems=f for problem in solver.problems has_converged = true if is_field_problem(problem) - has_converged = problem.assembly.u_norm_change/norm(problem.assembly.u) < eps + has_converged = problem.assembly.u_norm_change < eps if isapprox(norm(problem.assembly.u), 0.0) has_converged = true end @@ -395,8 +397,8 @@ function call(solver::Solver) # 2.3 update solution back to elements for problem in solver.problems - u, la = update_assembly!(problem, u, la) - update_elements!(problem, u, la) + u_new, la_new = update_assembly!(problem, u, la) + update_elements!(problem, u_new, la_new) end # 2.4 check convergence @@ -407,5 +409,7 @@ function call(solver::Solver) end # 3. did not converge - throw(NonlinearConvergenceError(solver)) + if solver.nonlinear_system_error_if_no_convergence + throw(NonlinearConvergenceError(solver)) + end end