fixed a lot of tests

src/JuliaFEM.jl

@@ -73,8 +73,16 @@ export Modal
 include("optics.jl")
 export find_intersection, calc_reflection, calc_normal
 
-### MORTAR STUFF ###
-include("mortar.jl")  # mortar projection
+### Mortar methods ###
+include("mortar.jl")
+export calculate_normals,
+       calculate_normals!,
+       project_from_slave_to_master,
+       project_from_master_to_slave,
+       Mortar
+
+### Contact mechanics ###
+#include("contact.jl")
 
 # rest of things
 include("utils.jl")

src/contact.jl

@@ -0,0 +1,130 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
+
+"""
+Currently two strategies exists:
+
+a) Remove inactive inequality constraints in element level. This is done in
+   assemble! if normal_condition is set to :Contact. For some reason this
+   leads to convergence issues.
+b) Remove inactive inequality constraints in assembly level. This is done in
+   posthook algorithm if inequality_constraints is set to true. This gives
+   more robust behavior.
+
+   Either use inequality_constraints=True OR :Contact + :Slip, but do not mix.
+
+   minimum_distance can be used to roughly skip integration of mortar
+   projections for elements that are "far enough" from each other. Increases
+   performance.
+
+"""
+type Mortar <: BoundaryProblem
+    formulation :: Symbol # :total, :incremental, :autodiff
+    dual_basis :: Bool
+    inequality_constraints :: Bool # Launch PDASS to solve inequality constraints
+    normal_condition :: Symbol # Tie or Contact
+    tangential_condition :: Symbol # Stick or Slip
+    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, false)
+end
+
+function get_unknown_field_name(::Type{Mortar})
+    return "reaction force"
+end
+
+function get_formulation_type(problem::Problem{Mortar})
+    return problem.properties.formulation
+end
+
+macro debug(msg)
+    haskey(ENV, "DEBUG") || return
+    return msg
+end
+
+function assemble!(problem::Problem{Mortar}, time::Real)
+    elements = get_elements(problem)
+    if length(elements) == 0
+        info("$(typeof(problem)) : forget to add elements?")
+        return
+    end
+    # returns 3 if eldim 2 (tri3, quad4, ...) for 3d problems etc.
+    eldim = size(elements[1], 1)+1
+    assemble!(problem, time, Val{eldim})
+end
+
+include("mortar_2d.jl")
+include("mortar_2d_autodiff.jl")
+include("mortar_3d.jl")
+include("mortar_3d_autodiff.jl")
+
+""" Remove inactive inequality constraints by using primal-dual active set strategy. """
+function boundary_assembly_posthook!(solver::Solver, problem::Problem{Mortar}, C1, C2, D, g)
+    problem.properties.inequality_constraints || return
+    info("PDASS: Starting primal-dual active set strategy to determine active constraints")
+    S = Set{Int64}()
+    for element in get_elements(problem)
+        haskey(element, "master elements") || continue
+        push!(S, get_connectivity(element)...)
+    end
+    S = sort(collect(S))
+    dim = get_unknown_field_dimension(problem)
+    ndofs = solver.ndofs
+    nnodes = round(Int, ndofs/dim)
+
+    c = reshape(full(problem.assembly.c, ndofs, 1), dim, nnodes)
+    A = find(c[1,:] .> 0)
+    A = intersect(A, S)
+    I = setdiff(S, A)
+
+    info("PDASS: contact nodes: $(sort(collect(S)))")
+    info("PDASS: active nodes: $(sort(collect(A)))")
+    info("PDASS: inactive nodes: $(sort(collect(I)))")
+
+    # remove any inactive nodes
+    for j in I
+        dofs = [dim*(j-1)+i for i=1:dim]
+        C1[dofs,:] = 0
+        C2[dofs,:] = 0
+        D[dofs,:] = 0
+        g[dofs,:] = 0
+    end
+
+    # handle tangential condition for active nodes
+    if problem.properties.tangential_condition == :Slip
+        for j in A
+            dofs = [dim*(j-1)+i for i=1:dim]
+            tangential_dofs = dofs[2:end]
+            D[tangential_dofs,dofs] = C2[tangential_dofs,dofs]
+            C2[tangential_dofs,:] = 0
+            g[tangential_dofs,:] = 0
+        end
+    end
+
+    return
+end
+
+function assemble_prehook!(problem::Problem{Mortar}, time::Real)
+    info("mortar assemble prehook at time $time")
+    slaves = Set{Element}()
+    for element in get_elements(problem)
+        haskey(element, "master elements") || continue
+        push!(slaves, element)
+    end
+    info("$(length(slaves)) slave elements")
+    length(slaves) != 0 || error("no slave elements found for problem (forget to add masters?).")
+    info("mortar: update normal-tangential system.")
+    calculate_normal_tangential_coordinates!(collect(slaves), time)
+    info("mortar assemble prehook done.")
+end

src/elasticity.jl

@@ -21,7 +21,6 @@ function get_unknown_field_name(problem::Problem{Elasticity})
 end
 
 function get_formulation_type(problem::Problem{Elasticity})
-    # we are solving residual and add increment to previous solution vector
     return :incremental
 end
 
@@ -57,21 +56,37 @@ function assemble{El<:Union{Tri3,Tri6,Quad4}}(problem::Problem{Elasticity}, elem
         N = element(ip, time)
         dN = element(ip, time, Val{:Grad})
 
-        # kinematics; calculate deformation gradient and strain
-        gradu = zeros(dim, dim)
-        if haskey(element, "displacement")
-            gradu += element("displacement", ip, time, Val{:Grad})
-        end
-        strain = zeros(dim , dim)
-        strain += 1/2*(gradu' + gradu)
-        F = eye(dim)
+        # kinematics
+
+        gradu = element("displacement", ip, time, Val{:Grad})
+        fill!(BL, 0.0)
+
         if props.finite_strain
-            F += gradu
-            strain += 1/2*gradu'*gradu
+            strain = 1/2*(gradu + gradu' + gradu'*gradu)
+            F = eye(dim) + gradu
+            for i=1:size(dN, 2)
+                BL[1, 2*(i-1)+1] += F[1,1]*dN[1,i]
+                BL[1, 2*(i-1)+2] += F[2,1]*dN[1,i]
+                BL[2, 2*(i-1)+1] += F[1,2]*dN[2,i]
+                BL[2, 2*(i-1)+2] += F[2,2]*dN[2,i]
+                BL[3, 2*(i-1)+1] += F[1,1]*dN[2,i] + F[1,2]*dN[1,i]
+                BL[3, 2*(i-1)+2] += F[2,1]*dN[2,i] + F[2,2]*dN[1,i]
+            end
+        else # linearized strain
+            strain = 1/2*(gradu + gradu')
+            F = eye(dim)
+            for i=1:size(dN, 2)
+                BL[1, 2*(i-1)+1] = dN[1,i]
+                BL[2, 2*(i-1)+2] = dN[2,i]
+                BL[3, 2*(i-1)+1] = dN[2,i]
+                BL[3, 2*(i-1)+2] = dN[1,i]
+            end
         end
 
-        # constitutive equations; material model (isotropic linear material here)
-        # get_material(problem, element, ...)
+        strain_vec = [strain[1,1]; strain[2,2]; strain[1,2]]
+        update!(ip, "strain", time => strain_vec)
+
+        # calculate stress
         E = element("youngs modulus", ip, time)
         nu = element("poissons ratio", ip, time)
         if props.formulation == :plane_stress
@@ -88,53 +103,48 @@ function assemble{El<:Union{Tri3,Tri6,Quad4}}(problem::Problem{Elasticity}, elem
             error("unknown plane formulation: $(props.formulation)")
         end
         # calculate stress
-        strain_vec = [strain[1,1]; strain[2,2]; 2.0*strain[1,2]]
-        stress_vec = D*strain_vec
-        stress = [stress_vec[1] stress_vec[3]; stress_vec[3] stress_vec[2]]
-        cauchy_stress = F'*stress*F/det(F)
-        cauchy_stress = [cauchy_stress[1,1]; cauchy_stress[2,2]; cauchy_stress[1,2]]
+        stress_vec = D * ([1.0, 1.0, 2.0] .* strain_vec)
+        update!(ip, "stress", time => stress_vec)
 
-        update!(ip, "strain", time => strain_vec)
-        update!(ip, "cauchy stress", time => cauchy_stress)
-        update!(ip, "pk2 stress", time => stress_vec)
+        Km += w*BL'*D*BL
 
-        # add contributions: material and geometric stiffness + internal forces
-        fill!(BL, 0.0)
-        for i=1:size(dN, 2)
-            BL[1, 2*(i-1)+1] = F[1,1]*dN[1,i]
-            BL[1, 2*(i-1)+2] = F[2,1]*dN[1,i]
-            BL[2, 2*(i-1)+1] = F[1,2]*dN[2,i]
-            BL[2, 2*(i-1)+2] = F[2,2]*dN[2,i]
-            BL[3, 2*(i-1)+1] = F[1,1]*dN[2,i] + F[1,2]*dN[1,i]
-            BL[3, 2*(i-1)+2] = F[2,1]*dN[2,i] + F[2,2]*dN[1,i]
-        end
-        fill!(BNL, 0.0)
-        for i=1:size(dN, 2)
-            BNL[1, 2*(i-1)+1] = dN[1,i]
-            BNL[2, 2*(i-1)+1] = dN[2,i]
-            BNL[3, 2*(i-1)+2] = dN[1,i]
-            BNL[4, 2*(i-1)+2] = dN[2,i]
-        end
-        S2 = zeros(2*dim, 2*dim)
-        S2[1,1] = stress_vec[1]
-        S2[2,2] = stress_vec[2]
-        S2[1,2] = S2[2,1] = stress_vec[3]
-        S2[3:4,3:4] = S2[1:2,1:2]
+        # stress = [stress_vec[1] stress_vec[3]; stress_vec[3] stress_vec[2]]
+        # cauchy_stress = F'*stress*F/det(F)
+        # cauchy_stress = [cauchy_stress[1,1]; cauchy_stress[2,2]; cauchy_stress[1,2]]
+        # update!(ip, "cauchy stress", time => cauchy_stress)
+        
+        # material stiffness end
 
-        Km += w*BL'*D*BL # material stiffness
-        if props.finite_strain # add geometric stiffness
+        if props.geometric_stiffness
+            # take geometric stiffness into account
+
+            fill!(BNL, 0.0)
+            for i=1:size(dN, 2)
+                BNL[1, 2*(i-1)+1] = dN[1,i]
+                BNL[2, 2*(i-1)+1] = dN[2,i]
+                BNL[3, 2*(i-1)+2] = dN[1,i]
+                BNL[4, 2*(i-1)+2] = dN[2,i]
+            end
+
+            S2 = zeros(2*dim, 2*dim)
+            S2[1,1] = stress_vec[1]
+            S2[2,2] = stress_vec[2]
+            S2[1,2] = S2[2,1] = stress_vec[3]
+            S2[3:4,3:4] = S2[1:2,1:2]
+            
             Kg += w*BNL'*S2*BNL # geometric stiffness
+
         end
 
-        if get_formulation_type(problem) == :incremental
-            f -= w*BL'*stress_vec # internal force
-        end
+        # rhs, internal and external load
+        
+        f -= w*BL'*stress_vec
 
-        # volume load
         if haskey(element, "displacement load")
             b = element("displacement load", ip, time)
             f += w*vec(N'*b)
         end
+
         for i=1:dim
             if haskey(element, "displacement load $i")
                 b = element("displacement load $i", ip, time)
@@ -409,7 +419,8 @@ function assemble{El<:Union{Tet4, Tet10, Hex8}}(problem::Problem{Elasticity}, el
 
         # material stiffness end
 
-        if props.geometric_stiffness # take geometric stiffness into account
+        if props.geometric_stiffness
+            # take geometric stiffness into account
 
             fill!(BNL, 0.0)
 

src/heat.jl

@@ -39,31 +39,32 @@ function get_unknown_field_type(problem::Problem{Heat})
 end
 
 function assemble!(assembly::Assembly, problem::Problem{Heat}, element::Element, time=0.0)
-
     gdofs = get_gdofs(problem, element)
-
+    field_name = get_unknown_field_name(problem)
+    nnodes = length(element)
+    K = zeros(nnodes, nnodes)
+    fq = zeros(nnodes)
     for ip in get_integration_points(element)
         detJ = element(ip, time, Val{:detJ})
         w = ip.weight*detJ
-
         N = element(ip, time)
-        if haskey(element, "density")
-            rho = element("density", ip, time)
-            add!(assembly.M, gdofs, gdofs, w*rho*N'*N)
-        end
-        if haskey(element, "temperature thermal conductivity")
+        if haskey(element, "$field_name thermal conductivity")
             dN = element(ip, time, Val{:Grad})
-            k = element("temperature thermal conductivity", ip, time)
-            add!(assembly.K, gdofs, gdofs, w*k*dN'*dN)
+            k = element("$field_name thermal conductivity", ip, time)
+            K += w*k*dN'*dN
         end
-        if haskey(element, "temperature load")
-            f = element("temperature load", ip, time)
-            add!(assembly.f, gdofs, w*N'*f)
+        if haskey(element, "$field_name load")
+            f = element("$field_name load", ip, time)
+            fq += w*N'*f
         end
-        if haskey(element, "temperature flux")
-            g = element("temperature flux", ip, time)
-            add!(assembly.f, gdofs, w*N'*g)
+        if haskey(element, "$field_name flux")
+            g = element("$field_name flux", ip, time)
+            fq += w*N'*g
         end
     end
+    T = vec(element[field_name](time))
+    fq -= K*T
+    add!(assembly.K, gdofs, gdofs, K)
+    add!(assembly.f, gdofs, fq)
 end
 

src/mortar.jl

@@ -1,43 +1,14 @@
 # This file is a part of JuliaFEM.
 # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
 
-"""
-Currently two strategies exists:
-
-a) Remove inactive inequality constraints in element level. This is done in
-   assemble! if normal_condition is set to :Contact. For some reason this
-   leads to convergence issues.
-b) Remove inactive inequality constraints in assembly level. This is done in
-   posthook algorithm if inequality_constraints is set to true. This gives
-   more robust behavior.
-
-   Either use inequality_constraints=True OR :Contact + :Slip, but do not mix.
-
-   minimum_distance can be used to roughly skip integration of mortar
-   projections for elements that are "far enough" from each other. Increases
-   performance.
-
-"""
 type Mortar <: BoundaryProblem
-    formulation :: Symbol # :total, :incremental, :autodiff
-    dual_basis :: Bool
-    inequality_constraints :: Bool # Launch PDASS to solve inequality constraints
-    normal_condition :: Symbol # Tie or Contact
-    tangential_condition :: Symbol # Stick or Slip
-    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
+    adjust :: Bool
+    tolerance :: Float64
 end
 
 function Mortar()
-    Mortar(:total, true, false, :Tie, :Stick, Inf, false, [], [], [], [], false, false, -1, false)
+    return Mortar(false, false, 0.0)
 end
 
 function get_unknown_field_name(::Type{Mortar})
@@ -45,86 +16,173 @@ function get_unknown_field_name(::Type{Mortar})
 end
 
 function get_formulation_type(problem::Problem{Mortar})
-    return problem.properties.formulation
+    return :incremental
 end
 
-macro debug(msg)
-    haskey(ENV, "DEBUG") || return
-    return msg
+typealias MortarElements2D Union{Seg2, Seg3}
+typealias MortarElements3D Union{Tri3, Tri6, Quad4}
+
+function newton(f, df, x; tol=1.0e-6, max_iterations=10)
+    for i=1:max_iterations
+        dx = -f(x)/df(x)
+        x += dx
+        if norm(dx) < tol
+            return x
+        end
+    end
+    error("Newton iteration did not converge in $max_iterations iterations")
 end
 
-function assemble!(problem::Problem{Mortar}, time::Real)
-    elements = get_elements(problem)
-    if length(elements) == 0
-        info("$(typeof(problem)) : forget to add elements?")
-        return
-    end
-    # returns 3 if eldim 2 (tri3, quad4, ...) for 3d problems etc.
-    eldim = size(elements[1], 1)+1
-    assemble!(problem, time, Val{eldim})
+function cross2(a, b)
+    cross([a; 0], [b; 0])[3]
 end
 
-include("mortar_2d.jl")
-include("mortar_2d_autodiff.jl")
-include("mortar_3d.jl")
-include("mortar_3d_autodiff.jl")
+function project_from_master_to_slave{E<:MortarElements2D}(slave_element::Element{E}, x2, time)
+    x1_ = slave_element["geometry"](time)
+    n1_ = slave_element["normal"](time)
+    x1(xi1) = vec(get_basis(slave_element, [xi1], time))*x1_
+    dx1(xi1) = vec(get_dbasis(slave_element, [xi1], time))*x1_
+    n1(xi1) = vec(get_basis(slave_element, [xi1], time))*n1_
+    dn1(xi1) = vec(get_dbasis(slave_element, [xi1], time))*n1_
+    R(xi1) = cross2(x1(xi1)-x2, n1(xi1))
+    dR(xi1) = cross2(dx1(xi1), n1(xi1)) + cross2(x1(xi1)-x2, dn1(xi1))
+    xi1 = newton(R, dR, 0.0)
+    return xi1
+end
 
-""" Remove inactive inequality constraints by using primal-dual active set strategy. """
-function boundary_assembly_posthook!(solver::Solver, problem::Problem{Mortar}, C1, C2, D, g)
-    problem.properties.inequality_constraints || return
-    info("PDASS: Starting primal-dual active set strategy to determine active constraints")
-    S = Set{Int64}()
-    for element in get_elements(problem)
-        haskey(element, "master elements") || continue
-        push!(S, get_connectivity(element)...)
-    end
-    S = sort(collect(S))
-    dim = get_unknown_field_dimension(problem)
-    ndofs = solver.ndofs
-    nnodes = round(Int, ndofs/dim)
+function project_from_slave_to_master{E<:MortarElements2D}(master_element::Element{E}, x1, n1, time)
+    x2_ = master_element["geometry"](time)
+    x2(xi2) = vec(get_basis(master_element, [xi2], time))*x2_
+    dx2(xi2) = vec(get_dbasis(master_element, [xi2], time))*x2_
+    cross2(a, b) = cross([a; 0], [b; 0])[3]
+    R(xi2) = cross2(x2(xi2)-x1, n1)
+    dR(xi2) = cross2(dx2(xi2), n1)
+    xi2 = newton(R, dR, 0.0)
+    return xi2
+end
 
-    c = reshape(full(problem.assembly.c, ndofs, 1), dim, nnodes)
-    A = find(c[1,:] .> 0)
-    A = intersect(A, S)
-    I = setdiff(S, A)
-
-    info("PDASS: contact nodes: $(sort(collect(S)))")
-    info("PDASS: active nodes: $(sort(collect(A)))")
-    info("PDASS: inactive nodes: $(sort(collect(I)))")
-
-    # remove any inactive nodes
-    for j in I
-        dofs = [dim*(j-1)+i for i=1:dim]
-        C1[dofs,:] = 0
-        C2[dofs,:] = 0
-        D[dofs,:] = 0
-        g[dofs,:] = 0
-    end
-
-    # handle tangential condition for active nodes
-    if problem.properties.tangential_condition == :Slip
-        for j in A
-            dofs = [dim*(j-1)+i for i=1:dim]
-            tangential_dofs = dofs[2:end]
-            D[tangential_dofs,dofs] = C2[tangential_dofs,dofs]
-            C2[tangential_dofs,:] = 0
-            g[tangential_dofs,:] = 0
+function calculate_normals(elements, time, rotate_normals=false)
+    tangents = Dict{Int64, Vector{Float64}}()
+    for element in elements
+        conn = get_connectivity(element)
+        X1 = element("geometry", time)
+        dN = get_dbasis(element, [0.0], time)
+        tangent = vec(sum([kron(dN[:,i], X1[i]') for i=1:length(X1)]))
+        for nid in conn
+            if haskey(tangents, nid)
+                tangents[nid] += tangent
+            else
+                tangents[nid] = tangent
+            end
         end
     end
 
-    return
+    Q = [0.0 -1.0; 1.0 0.0]
+    normals = Dict{Int64, Vector{Float64}}()
+    S = sort(collect(keys(tangents)))
+    for j in S
+        tangents[j] /= norm(tangents[j])
+        normals[j] = Q*tangents[j]
+    end
+
+    if rotate_normals
+        for j in S
+            normals[j] = -normals[j]
+        end
+    end
+
+    return normals, tangents
 end
 
-function assemble_prehook!(problem::Problem{Mortar}, time::Real)
-    info("mortar assemble prehook at time $time")
-    slaves = Set{Element}()
-    for element in get_elements(problem)
-        haskey(element, "master elements") || continue
-        push!(slaves, element)
+function calculate_normals!(elements, time, rotate_normals=false)
+    normals, tangents = calculate_normals(elements, time, rotate_normals)
+    for element in elements
+        conn = get_connectivity(element)
+        update!(element, "normal", time => [normals[j] for j in conn])
+        update!(element, "tangent", time => [tangents[j] for j in conn])
     end
-    info("$(length(slaves)) slave elements")
-    length(slaves) != 0 || error("no slave elements found for problem (forget to add masters?).")
-    info("mortar: update normal-tangential system.")
-    calculate_normal_tangential_coordinates!(collect(slaves), time)
-    info("mortar assemble prehook done.")
+end
+
+function assemble!(problem::Problem{Mortar}, time::Real)
+
+    props = problem.properties
+    field_dim = get_unknown_field_dimension(problem)
+    field_name = get_parent_field_name(problem)
+    slave_elements = filter(el -> haskey(el, "master elements"), get_elements(problem))
+
+    # 1. calculate nodal normals and tangents for slave element nodes j ∈ S
+    normals, tangents = calculate_normals(slave_elements, time, props.rotate_normals)
+    update!(slave_elements, "normal", normals)
+    update!(slave_elements, "tangent", tangents)
+    S = sort(collect(keys(normals)))
+
+    # 2. loop all slave elements
+    for slave_element in slave_elements
+        haskey(slave_element, "master elements") || continue
+
+        slave_element_nodes = get_connectivity(slave_element)
+        nsl = length(slave_element)
+        X1 = slave_element["geometry"](time)
+        n1 = Field([normals[j] for j in slave_element_nodes])
+
+        # 3. loop all master elements
+        for master_element in slave_element["master elements"](time)
+
+            master_element_nodes = get_connectivity(master_element)
+            nm = length(master_element)
+            X2 = master_element["geometry"](time)
+
+            # 3.1 calculate segmentation
+            xi1a = project_from_master_to_slave(slave_element, X2[1], time)
+            xi1b = project_from_master_to_slave(slave_element, X2[end], 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.3. loop integration points of one integration segment and calculate
+            # local mortar matrices
+            De = zeros(nsl, nsl)
+            Me = zeros(nsl, nm)
+            ge = zeros(nsl)
+            for ip in get_integration_points(slave_element, 2)
+                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))
+                # 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*N1*N1'
+                Me += w*N1*N2'
+                if props.adjust
+                    g = X_s-X_m
+                    if g < props.tol
+                        ge += w*g
+                    end
+                end
+            end
+
+            # add contribution to contact virtual work
+            sdofs = get_gdofs(problem, slave_element)
+            mdofs = get_gdofs(problem, master_element)
+            
+            for i=1:field_dim
+                lsdofs = sdofs[i:field_dim:end]
+                lmdofs = mdofs[i:field_dim:end]
+                add!(problem.assembly.C1, lsdofs, lsdofs, De)
+                add!(problem.assembly.C1, lsdofs, lmdofs, -Me)
+                add!(problem.assembly.C2, lsdofs, lsdofs, De)
+                add!(problem.assembly.C2, lsdofs, lmdofs, -Me)
+                add!(problem.assembly.g, lsdofs, ge)
+            end
+
+
+        end # master elements done
+
+    end # slave elements done, contact virtual work ready
+
 end

src/problems.jl

@@ -111,11 +111,11 @@ function Problem{P<:BoundaryProblem}(::Type{P}, main_problem::Problem, elements=
 end
 
 function get_formulation_type{P<:FieldProblem}(problem::Problem{P})
-    return :total
+    return :incremental
 end
 
 function get_formulation_type{P<:BoundaryProblem}(problem::Problem{P})
-    return :total
+    return :incremental
 end
 
 function get_assembly(problem)

src/solvers.jl

@@ -256,7 +256,7 @@ function solve_linear_system(solver::Solver, ::Type{Val{:DirectLinearSolver}})
     # assemble boundary problems
     Kb, C1, C2, D, fb, g = get_boundary_assembly(solver)
 
-    K = K + Kb
+    K = K + Kb + Kg
     f = f + fb
     K = 1/2*(K + K')
     u = zeros(solver.ndofs)
@@ -279,7 +279,7 @@ function solve_linear_system(solver::Solver, ::Type{Val{:DirectLinearSolver}})
     end
 
     # solver interior
-    CF = cholfact(K[interior_dofs, interior_dofs])
+    CF = ldltfact(K[interior_dofs, interior_dofs])
     Kib = K[interior_dofs, boundary_dofs]
     Kbb = K[boundary_dofs, boundary_dofs]
     fi = f[interior_dofs]
@@ -337,6 +337,19 @@ function Base.showerror(io::IO, exception::NonlinearConvergenceError)
     print(io, "nonlinear iteration did not converge in $max_iters iterations!")
 end
 
+function assemble!(solver::Solver; force_assembly=true)
+    info("Assembling problems ...")
+    tic()
+    for problem in solver.problems
+        if force_assembly # force reassembly
+            problem.assembly.changed = true
+        end
+        assemble!(problem, solver.time)
+    end
+    t1 = round(toq(), 2)
+    info("Assembled in $t1 seconds.")
+end
+
 """ Default solver for quasistatic nonlinear problems. """
 function call(solver::Solver{Nonlinear})
 
@@ -352,14 +365,7 @@ function call(solver::Solver{Nonlinear})
         info("Starting nonlinear iteration #$(properties.iteration)")
 
         # 2.1 update linearized assemblies (if needed)
-        info("Assembling problems ...")
-        tic()
-        for problem in solver.problems
-            problem.assembly.changed = true  # force reassembly
-            assemble!(problem, solver.time)
-        end
-        t1 = round(toq(), 2)
-        info("Assembled in $t1 seconds.")
+        assemble!(solver)
 
         # 2.2 call solver for linearized system (default: direct lu factorization)
         info("Solve linear system ...")

src/sparse.jl

@@ -155,6 +155,10 @@ function get_nonzero_rows(A::SparseMatrixCOO)
     return get_nonzero_rows(sparse(A))
 end
 
+function get_nonzero_rows(A::Matrix)
+    return get_nonzero_rows(sparse(A))
+end
+
 function size(A::SparseMatrixCOO)
     return maximum(A.I), maximum(A.J)
 end