From 64becfb29dccccbb64dd93112cc71e71a9ea951c Mon Sep 17 00:00:00 2001
From: Jukka Aho <jukka.aho@kapsi.fi>
Date: Mon, 3 Oct 2016 23:12:00 +0300
Subject: [PATCH] 3d finite sliding contact

---
 src/problems_contact.jl   | 22 +++++++++++++++++++++-
 src/problems_mortar_2d.jl |  3 ++-
 src/problems_mortar_3d.jl | 20 ++++++++++++--------
 3 files changed, 35 insertions(+), 10 deletions(-)

diff --git a/src/problems_contact.jl b/src/problems_contact.jl
index f15abef..483bfeb 100644
--- a/src/problems_contact.jl
+++ b/src/problems_contact.jl
@@ -17,9 +17,14 @@ type Contact <: BoundaryProblem
     friction :: Bool
     dual_basis :: Bool
     use_forwarddiff :: Bool
+    forwarddiff_assemble_in_pieces :: Bool
     minimum_active_set_size :: Int
     distval :: Float64
     remove_from_set :: Bool # allow removal of non-potential contact pairs
+    allow_quads :: Bool
+    remove_nodes :: Vector{Int}
+    always_in_contact :: Bool
+    update_contact_pairing :: Bool
     store_fields :: Vector{AbstractString}
 end
 
@@ -27,7 +32,22 @@ function Contact()
     default_fields = ["element area", "contact area", "weighted gap",
         "contact pressure", "active nodes", "inactive nodes", "stick nodes",
         "slip nodes", "complementarity condition", "contact error"]
-    return Contact(-1, false, false, false, true, false, 0, 5.0, false, default_fields)
+    return Contact(
+        -1,      # dimension
+        false,   # rotate_normals
+        false,   # finite_sliding
+        false,   # friciton
+        true,    # dual basis
+        false,   # use forwarddiff
+        false,   # when using forwarddiff, assemble interface in pieces
+        0,       # minimum active set size
+        5.0,     # distance value for contact detection
+        false,   # allow removal of non-potential contact pairs
+        false,   # allow quadrangles in contact discretization
+        [],      # remove these nodes always from set
+        false,   # mainly for debugging, do not remove inactive nodes
+        true,    # update contact pairing on each loop
+        default_fields)
 end
 
 function get_unknown_field_name(problem::Problem{Contact})
diff --git a/src/problems_mortar_2d.jl b/src/problems_mortar_2d.jl
index 7ddf890..19c830a 100644
--- a/src/problems_mortar_2d.jl
+++ b/src/problems_mortar_2d.jl
@@ -19,7 +19,8 @@ function cross2(a, b)
 end
 
 function get_slave_elements(problem::Problem)
-    filter(el -> haskey(el, "master elements"), get_elements(problem))
+    cond(el) = haskey(el, "master elements") || haskey(el, "potential master elements")
+    return filter(cond, get_elements(problem))
 end
 
 function project_from_master_to_slave{E<:MortarElements2D}(slave_element::Element{E}, x2, time)
diff --git a/src/problems_mortar_3d.jl b/src/problems_mortar_3d.jl
index 8ae3699..1fb5822 100644
--- a/src/problems_mortar_3d.jl
+++ b/src/problems_mortar_3d.jl
@@ -23,7 +23,7 @@ function inv3(P::Matrix)
     return 1/(a*A + b*B + c*C)*[A B C; D E F; G H I]
 end
 
-function vertex_inside_polygon(q, P; atol=1.0e-6)
+function vertex_inside_polygon(q, P; atol=1.0e-3)
     N = length(P)
     angle = 0.0
     for i=1:N
@@ -56,14 +56,18 @@ function calculate_centroid(P)
     return C
 end
 
-function get_cells(P, C)
+function get_cells(P, C; allow_quads=false)
     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]
+    # trivial cases, polygon already triangle / quadrangle
+    if N == 3
+        return Vector[P]
+    end
+    if N == 4 && allow_quads
+        return Vector[P]
+    end
     #V = sum([cross(P[i], P[mod(i,N)+1]) for i=1:N])
     #A = 1/2*abs(dot(n, V))
     #info("A = $A")
@@ -87,7 +91,7 @@ function get_cells(P, C)
 end
 
 """ Test does P contain q. """
-function contains{T}(P::Vector{T}, q::T; check_is_close=true, rtol=1.0e-5)
+function contains{T}(P::Vector{T}, q::T; check_is_close=true, rtol=1.0e-4)
     if q in P
         return true
     end
@@ -105,7 +109,7 @@ function get_polygon_clip(xs, xm, n; debug=false)
     # objective: search does line xm1 - xm2 clip xs
     nm = length(xm)
     ns = length(xs)
-    P = Vector{Number}[]
+    P = Vector[]
 
     # 1. test is master point inside slave, if yes, add to clip
     for i=1:nm
@@ -152,7 +156,7 @@ function get_polygon_clip(xs, xm, n; debug=false)
 end
 
 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)
+    element::Element{E}, x::DVTI, time::Real; max_iterations::Int=10, iter_tol::Float64=1.0e-6)
     basis(xi) = get_basis(element, xi, time)
     dbasis(xi) = get_dbasis(element, xi, time)
     f(theta) = basis(theta[1:2])*x - theta[3]*n0 - p