Incremental formulation works now, also finite sliding without proper linearization (very slow convergence)

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

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
 

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))

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
-

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
-

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