.. _contact:

#######
Contact
#######

\section{Contact virtual work}

As a starting point for the treatment of contact, its contribution to the virtual work expression can be stated as:

.. math::
   :label: contact:eq:01

   \int_{S^3} \left( -t_N \delta g_N + t_{T_\alpha} \delta g_T^{\alpha} \right) \text{d} a ,


where :math:`S^3` is the common surface between two continua, :math:`t_N` is the contact normal traction (positive in compression)
:math:`t_{T_\alpha}` is the contact tangential traction in one of two local (tangent plane) directions :math:`\alpha`, and $\delta
g_N:math:` and `\delta g_T^{\alpha}:math:` are the directional derivatives of the contact normal gap  `g_N$ and tangential slip
:math:`g_T^{\alpha}` in the direction of :math:`\dot{\varphi}`, i.e.:

.. math::
   :label: contact:eq:02

   \begin{aligned}
   \delta g_N &:= \left.\frac{\textrm{d}}{\textrm{d}\beta}\right|_{\beta=0} \left[ g_N \left(\varphi +\beta \dot{\varphi}\right) \right] \\,
   \delta g_T^{\alpha} &:= \left.\frac{\textrm{d}}{\textrm{d}\beta}\right|_{\beta=0} \left[ g_T^{\alpha} \left(\varphi +\beta \dot{\varphi}\right) \right] .
   \end{aligned}

In :eq:`contact:eq:01`, the deformation is subject to the following constraints, referred to as the **Kuhn-Tucker conditions**. The Kuhn-Tucker conditions are a set of constraints to be considered representative of the mechanical contact problem in continuum mechanics, and can be written as:

.. math::
   :label: contact:eq:03

    \begin{aligned}
    \delta g_N \geq 0 & \text{  (a) impenetrability constraint},\\
    t_N \geq 0            & \text{  (b) no adhesion condition},\\
    t_N g_n \geq 0      & \text{  (c) complementary condition}, \\
    t_N \dot{g}_N \geq 0 & \text{  (d) persistency condition},
    \end{aligned}

for frictionless response and

.. math::
   :label: contact:eq:04

    \begin{aligned}
    \Phi := \| t_T \| - \mu t_N \leq 0 & \text{  (a) slip function}, \\
    L_v g_T -\varsigma \frac{t_T}{\|t_T\|} = 0 & \text{  (b) slip rule},\\
    \varsigma \geq 0 & \text{  (c) consistency parameter},\\
    \Phi \varsigma = 0 & \text{  (d) complementary condition},
    \end{aligned}


for the prescription of a Coulomb friction (where :math:`\mu` is the friction coefficient).
In :eq:`contact:eq:03`, the gap :math:`g_N` is defined with respect to all material points :math:`\mathbf{Y} \in S^3` as:

.. math::
   :label: contact:eq:05

   g_N ( \mathbf{X},t) = \min_{\mathbf{Y} \in S^3} \| \varphi ( \mathbf{X} ) - \varphi ( \mathbf{Y} ) \| *sign* (g_N),

where

.. math::
   :label: contact:eq:06

   \textrm{sign} (g_N) = \begin{cases}
    -1 & \text{if } \varphi (\mathbf{X}) \text{ lies on the interior of the contacted body}, \\
    1 & \text{otherwise} .
    \end{cases}

The tangential gap rate :math:`L_v g_T` in :eq:`contact:eq:04` is defined as follows:

.. math::
   :label: contact:eq:07

    L_v g_T = \left( \varphi \left( \mathbf{X} \right) - \varphi \left( \bar{\mathbf{Y}} (\mathbf{X}) \right) \right) \cdot
    \left( \mathbf{p}^{\alpha} \otimes \mathbf{p}_{\alpha} \right),


where :math:`\mathbf{p}^{\alpha}` and :math:`\mathbf{p}_{\alpha}` are base vectors associated with any appropriate surface
coordinate system used to describe :math:`S^3`, with these base vectors being evaluated at the *current* contact point
:math:`(\bar{\mathbf{Y}}(\mathbf{X}))` that satisfies the minimization of :eq:`contact:eq:03`. Use of the notation
:math:`L_v` is meant to imply a Lie derivative, which can be understood to be the time derivative of an object as viewed from
an embedded reference frame, in this case the convected frame :math:`\mathbf{p}_{\alpha}` frame, that moves along with the
point.

****************************************
Discretized forms of contact constraints
****************************************

The question is then, how to represent these conditions in discretized form suitable for FE solution methods. A simple example, shown in :numref:`contact-fig-constraint_types`, serves to demonstrate the concern that this question embodies. Two discretizations for the interface are evident and, as this simple example indicates, leads to an ambiguous definition of the interface.  

.. _contact-fig-constraint_types:

.. figure:: ../_static/figures/17-1.png
   :align: center
   :scale: 75 %

   Concerns in constraints choices for contact problems.

A historical treatment of contact has focused on applying the Kuhn-Tucker condition directly to the discretized form,
leading to what we will be referring to here as a node-face treatment of contact, or **node-face contact**. As we
will review here in this chapter, this treatment of contact is relatively straightforward from a conceptual standpoint,
however it does have several issues - even to the point of the overall approach being pathological in some applications.  

Alternatively, more recent investigations have focused on addressing these issues, leading to what we will be referring
to here as a face-face treatment of contact, or **face-face contact**. These methods consider the weak form more
directly, thus leading to a variationally consistent approach (e.g., mortar methods are an example of this approach and
are, at the moment, prevalent in the literature).

Node-Face contact
=================

For node-face contact the Kuhn-Tucker conditions are assumed to apply to one side of the contacting surfaces.Thus the gap :math:`g_N` is defined with respect to all **nodal points** :math:`\mathbf{Y}_I` as:

.. math::
   :label: contact:eq:08

    g_N ( \mathbf{X}_I,t) = \min_{\mathbf{Y}_I \in S^3} \| \varphi ( \mathbf{X}_I ) - \varphi ( \mathbf{Y} ) \|
    *sign* (g_N),


where :math:`I` refers to a nodal point on one side of the interface, whose coordinates are :math:`\mathbf{X}_I,t` at time :math:`t` of interest.
The right-hand side of :eq:`contact:eq:08` is the discrete form of :eq:`contact:eq:05` but is more commonly called the **closest-point projection**, which will be discussed in some detail in :numref:`contact-contact_kinematics`. 

As mentioned, there are issues associated with node-face constraints.  They stem from the application of contact
constraints directly to the discretized problem.  As shown in :numref:`contact-fig-constraint_types`, the potential
to over constrain the interface is avoided by applying the impenetrability constraint only at selected points along the
interface.  In the Solid Mechanics module these points coincide with the nodes, as it makes it convenient to obtain
contact results (normal and tangential tractions, stick/slip results, etc.) and interpret them in post-processing.

However, this approach does not truly alleviate over-constraining.  This is
easily demonstrated with an enlightening example (we will make use of this
example for the discussion of node-face contact and face-face contact, so
making a proper introduction is worthwhile).
:numref:`contact-fig-beam_pure_bending`, shows a beam bending problem that
is being modeled with continuum elements, in this case hex8 elements.  The
beam is cantilevered at its left end appropriately, i.e., fixed at the neutral
axis and constrained from motion only in the x-direction elsewhere. 

.. _contact-fig-beam_pure_bending:

.. figure:: ../_static/figures/17-2.png
   :align: center
   :scale: 75 %

   A continuum beam subjected to pure bending.

The analytic solution to this problem is one where the neutral axis should take the displacement corresponding to an arc
of a circle. When the moment is prescribed to be :math:`M^*` the beam should deform into a perfect circle.

When the beam is meshed with either an all coarse mesh (4 elements through its thickness) or an all fine mesh (16
elements through its thickness), the Finite Element results appear to be quite acceptable, producing the pure bending
solution.

However, lets now combine coarse and fine discretizations to solve the problem. In this case a mesh tying constraint is
required to obtain the solution, *which is seen to be fundamentally a contact problem with adhesion and infinite frictional capacity*.  The combining of coarse and fine discretizations can be done is a couple of canonical ways, as
shown in :numref:`contact-fig-beam_pure_bending2`; one where the interface between the discretizations is vertical
and the other where is along the neutral axis.

.. _contact-fig-beam_pure_bending2:

.. figure:: ../_static/figures/17-3.png
   :align: center
   :scale: 75 %

   A continuum beam that includes mesh tying subjected to pure bending.

In both cases, we apply the standard rule of thumb: given the same material on both sides of the interface, apply the
contact constraints on the finer discretization. Subjecting the beam to the prescribed moment reveals at once the issue:
kinematically enforcing a zero gap condition at each node is exactly correct in one case, where the interface is through
the depth of the beam, and severely over constraining in the other, where the interface is along the neutral axis. As
:numref:`contact-fig-beam_pure_bending_MPC` shows, the over constraint can be severe and may produce spurious stress
distributions in the fine mesh, particularly near the neutral axis.

.. _contact-fig-beam_pure_bending_MPC:

.. figure:: ../_static/figures/17-4.png
   :align: center
   :scale: 75 %

   Results for a continuum beam that includes mesh tying subjected to pure bending.

.. %\subsection{Face-Face contact}
.. %Face-face contact considers the weak form more directly. Examination of :eq:`contact:equation:01` yields an interpretation that :math:`t_N` as a Lagrange multiplier is meant to (weakly) enforce the kinematic quantity :math:` \delta g_N` to be zero.
.. %Along these lines, mortar methods have been developed and studied at length in the literature, see e.g., [**refs**].

.. %... thus leading to a variationally consistent approach.

**************
Contact Search
**************

The contact search algorithm is a logical component of the overall contact capability. Much of the reason to consider
search as a separate component is due to a need to revisit and replace algorithms as they demonstrate better
performance. As problem sizes grow there is an increasing computational cost of this aspect of computational solid
mechanics, particularly on distributed memory (parallel) computers.

As a way of introducing the concepts inherent in contact search algorithms, we recognize the similarity of contact
search to many other other problems in the simulation domain (e.g., the video gaming industry). In this more abstract
sense, a significant part of the contact search algorithm is a proximity determination of one object with respect to
another. Collision detection, as it is also referred to, is computationally intensive but also studied thoroughly to
obtain the best performance possible. Thus the contact algorithm in Sierra/SolidMechanics is comprised of the more general
proximity search followed by the much more specific detailed detection of contact in the context of a discrete Finite
Element method.

Proximity search algorithms
===========================

Although various proximity search algorithms have been developed over the years, those that have been used in Sierra/SolidMechanics are discussed.  Proximity search algorithms deal with bounding boxes. Construction of bounding boxes are
straightforwardly computed as the vector of min/max coordinates of the volume swept by the predicted motion of either
nodes or faces over a time step.  Specifically, over the time step :math:`( t \rightarrow t+\Delta t )`, the axis-aligned
bounding box for node :math:`I` is computed as follows: 

.. math::
   :label: contact:eq:09

    aabb(N_I) = \min( \varphi(\mathbf{X}_I , t) , \varphi(\mathbf{X}_I , t+\Delta t) ) \text{ , } \max( \varphi(\mathbf{X}_I
    , t) , \varphi(\mathbf{X}_I , t+\Delta t) ).


The resulting data of an axis aligned bounding box calculation is the absolute minimum amount of data (min and max
:math:`x,y,z` coordinates) representing a contact entity. The example expressed in :eq:`contact:eq:09` is for a
node; with the extension to a face being straightforwardly computed as:

.. math::
   :label: contact:eq:10

    aabb(F_M) = \min( aabb(N_J), \text{ }J=\text{nodes of face } M )\text{ , } \max( aabb(N_J), \text{ }J=\text{nodes of
    face } M ).


This allows a simplification or reduction of the contact problem to the more
general proximity detection based solely on axis-aligned bounding boxes of
contact entities.  Many of the structural modeling capabilities within
Sierra/SM are converted to these primitives (i.e., nodes and faces). Beams and
shells, for example, do not explicitly represent volume, but a volume is
inferred with ancillary data such as thickness.  These elements are converted
to contact primitives by lofting the finite element geometry explicitly to
volumetric discretizations.

.. %**Bucket search** 

.. %This is an algorithm, reference[],  that uniformly divides the volume of interest into volumes (buckets), places all contact entities into buckets as a coarse-grain search, and then does an :math:`n^2` search among all objects in a bucket.
.. %ADD REFERENCE.

.. %**Binary space search** 

.. %This is an algorithm, reference [], that is a serial search algorithm.



.. %** KDTree **

.. %This is the search algorithm that is used in SM.


Parallel search algorithms
==========================

The strategy for computational simulation of contact on distributed memory architectures (parallel computing) is to
decompose the contact problem in a distributed manner among the compute nodes. Typically known as domain decomposition
methods, there are several that are directly applicable to the contact problem. Inertial decomposition and recursive
coordinate bisection (RCB) decomposition are two geometric-based algorithms that are examples.  This geometry-based
decomposition approach is depicted graphically in :numref:`contact-fig-search1`.

.. _contact-fig-search1:

.. figure:: ../_static/figures/17-5.png
   :align: center
   :scale: 75 %

   Simple illustration of the domain decomposition for contact problems.

The proximity search algorithm thus performs with parallel scalability and serial efficiently on each processor. 

.. _contact-contact_kinematics:

Contact kinematics
==================

Recall that the output of the proximity search is a collection of bounding box pairs whose volumes overlap. The contact
entities associated with the bounding boxes are then considered for contact in what we call a detailed search. Detailed
searching is a term applied to computing the contact kinematic quantities associated with either the node-face or
face-face algorithm. 

**Closest point projection for node-face contact** 

The underpinning of the node-face contact approach is the choice of a set of points at which to apply the Kuhn-Tucker
conditions. Once this choice is made, it follows that a contact point must be determined for each node. The closest
point projection is the name given to this calculation, which is simply the point on the opposing surface that minimizes
the gap, i.e.,

.. math::
   :label: contact:eq:11

    g_N ( \mathbf{X_I},t) = \min_{\mathbf{Y} \in \{ *facets on side A surface* \} } \| \varphi ( \mathbf{X_I} ) -
    \varphi ( \mathbf{Y} ) \| *sign* (g_N).


This can be seen simply as the discrete form of the gap function
expressed in :eq:`contact:eq:03`.  Immediately with the
discrete surface *not* possessing :math:`C_1` continuity, the closest
point projection is no longer unique. At and around edges and corners
are the regions on the side A surface where the issues arise. :numref:`contact-fig-cpp_not_unique` depicts a simple illustration of the non-uniqueness encountered.

.. _contact-fig-cpp_not_unique:

.. figure:: ../_static/figures/17-6.png
   :align: center
   :scale: 75 %

   Simple illustration of the non-uniqueness in the closest point projection.

**Minimum volume overlap for face-face contact** 

The face-face contact approach seeks to avoid these issues by considering the volume overlap between the discrete sides of the interface. 


