# Contact Boundary Conditions

## Formulation
Rigid-body contact boundary conditions are implemented in Ratel using a Nitsche method approach as described in {cite}`mlika2018nitsche`, which we describe in detail for completeness.
Currently, only a rigid plane (platen) is supported.

### Large-Deformation Frictional Contact

Given a boundary region $\Gamma_c$, the gap between the rigid shape $\mathcal{B}$ and a deformed configuration point $\bm{x}\in\Gamma_c$ is given by the distance between $\bm{x}$ and the nearest point on the shape $\bm{y}\in\partial\mathcal{B}$.
That is,

$$
    \bm{y} = \arg\min_{\bm y\in\partial\mathcal{B}} \| x - y \|_2,
$$

and the signed gap function is given by

$$
   g(\bm{x}) = \bm{n}_y \cdot \left( \bm{x} - \bm{y} \right).
$$

We can formulate the gap in terms of the displacement field $\bm u$ as

$$
    g(\bm u) = \bm{n}_y \cdot(\bm{u} + \bm X - \bm c) = \bm{n}_y \cdot u + \bm{n}_y (\bm X - \bm c) = \bm{n}_y \cdot u + g_0,
$$

where $g_0$ is the gap in the initial configuration.

Denote the contact traction by $\hat{\bm{\sigma}}_N(\bm{u})$, which is computed as the first Piola-Kirchoff traction vector $\bm{P} \cdot \hat{\bm N}$, where $\hat{\bm N}$ is the normal at the material point $\bm{X} = \phi^{-1}(\bm{x})$.
Note that regardless of whether the volume operator is computed using the initial or current configurations, the contact traction is always computed in the initial configuration to simplify the linearization.
The contact traction is then decomposed into normal and tangential components, relative to the normal $\bm{n}_y$, as

$$
    \hat{\bm{\sigma}}_{N}(\bm{u}) = -\hat{\sigma}_n (u) \bm{n}_y + \bm{T}_{\bm{n}_y}\hat{\bm{\sigma}}_{N}(\bm{u}) = -\hat{\sigma}_n (u) \bm{n}_y + \hat{\bm{\sigma}}_t,
$$

where $\bm{T}_{\bm{v}} = \bm{I} - \bm{v} \otimes \bm{v}$ is the projection operator onto the tangent plane defined by the normal vector $\bm{v}$.
The normal component $\hat{\sigma}_n$ is the contact pressure.

Using these quantities, normal contact requires the following:

$$
\begin{aligned}
    g(\bm u) &\geq 0\\
    \hat{\sigma}_n (\bm u) &\leq 0\\
    \hat{\sigma}_n g(\bm u) &= 0.
\end{aligned}
$$ (normal-contact)

That is, if the point is in contact with the shape, then the contact pressure must be negative; otherwise, the contact pressure must be zero.

For Coulomb friction, we also require:

$$
\begin{cases}
    \|\hat{\bm{\sigma}}_t\| \leq -\mathscr{F} \hat{\sigma}_n(\bm u) & \text{if } \dot{\bm{u}} = 0,\\
    \hat{\bm{\sigma}}_t = \mathscr{F} \hat{\sigma}_n(\bm u) \frac{\dot{\bm{u}}}{\|\dot{\bm{u}}\|} & \text{otherwise}.
\end{cases}
$$ (coulomb-friction-contact)

Adopting the notation of {cite}`mlika2018nitsche`, we define the projection operator $P_{B(\bm{n}, \delta)}(\bm{v})$ which first projects the vector $\bm{v}$ onto the tangent plane defined by the normal vector $\bm{n}$, then projects the result onto a ball of radius $\delta$.
Specifically, that is,

$$
    P_{B(\bm{n}, \delta)}(\bm{v}) = \begin{cases}
        \bm{T}_{\bm{n}} \bm{v} & \text{if } \|\bm{T}_{\bm{n}} \bm{v}\| \leq \delta, \\
        \delta \frac{\bm{T}_{\bm{n}} \bm{v}}{\|\bm{T}_{\bm{n}} \bm{v}\|} & \text{otherwise}.
    \end{cases}
$$

Additional friction models are also supported, see [Friction Models](friction).

### Nitsche Method
The Nitsche method reformulates the conditions {eq}`normal-contact` as

$$
    \hat{\sigma}_n (\bm{u}) = \left[ \hat{\sigma}_n(\bm u) + \gamma g(\bm u) \right]_{\mathbb{R}^-},
$$ (nitsche-normal-force)

where $[v]_{\mathbb{R}^-} = \frac{1}{2}(v-|v|)$ is the projection onto the negative real numbers and $\gamma > 0$ is the Nitsche method parameter.
The Nitsche method parameter $\gamma$ serves a similar role to a penalty method, and larger values of $\gamma$ result in more robust non-penetration.
In general, $\gamma$ need only be a positive function defined on $\Gamma_c$.

The friction conditions {eq}`coulomb-friction-contact` are similarly reformulated as

$$
    \hat{\bm{\sigma}}_t = P_{B\left(\bm{n}_y, -\mathscr{F\left[ \hat{\sigma}_n(\bm u) + \gamma g(\bm u) \right]_{\mathbb{R}^-}}\right)}(\hat{\bm{\sigma}}_N - \gamma \dot{\bm{u}}),
$$

The resulting formulation for frictional contact is then

$$
\hat{\bm{\sigma}}_N = -\left[ -\hat{\bm{\sigma}}_N \cdot \bm{n}_y + \gamma g \right]_{\mathbb{R}^-}\bm{n}_y + \left[\hat{\bm{\sigma}}_N - \gamma \dot{\bm{u}}\right]_{-\mathscr{F}\left[ -\hat{\bm{\sigma}}_N\cdot\bm{n}_y + \gamma g \right]_{\mathbb{R}^-}}.
$$

### Weak Form

For simplicity, we only describe the boundary part of the weak form, which replaces the final (boundary) term in, e.g. {eq}`hyperelastic-weak-form-initial`.

There are several variants of the Nitsche method formulation depending on an auxiliary parameter $\theta\in\mathbb{R}$.
Typically, $\theta=-1,0,1$ are the values considered.
Currently, Ratel implements boundary conditions with $\theta=0$, which results in the simplest Jacobian.
The weak form of the boundary conditions on $\Gamma_c$ can be written as:

$$
    \int_{\Gamma_c} \bm{v}\cdot\left(\left[ -\hat{\bm{\sigma}}_N \cdot \bm{n}_y + \gamma g \right]_{\mathbb{R}^-}\bm{n}_y\right) \ dS - \int_{\Gamma_c} \bm{v}\cdot P_{B\left(\bm{n}_y, -\mathscr{F}\left[ -\hat{\bm{\sigma}}_N\cdot\bm{n}_y + \gamma g \right]_{\mathbb{R}^-}\right)}(\hat{\bm{\sigma}}_N - \gamma \dot{\bm{u}}) \ dS
$$ (contact-weak-form)

(contact-newton-linearization)=

### Newton Linearization

In order to solve using Newton-Krylov solvers, we require the linearization of {eq}`contact-weak-form`.

The linearization of the frictionless case is straightforward to compute.
Let $H$ be the Heaviside function, defined as

$$
    H(x) = \begin{cases}
        0, & x < 0;\\
        1, & x \geq 0.
    \end{cases}
$$

Then, we can compute

$$
    \mathrm{d}\left(\left[ \hat{\sigma}_n(\bm u) + \gamma g(\bm u) \right]_{\mathbb{R}^-}\right) = H\left(-\hat{\sigma}_n(\bm u) - \gamma g(\bm u)\right)
    \left( -\mathrm{d}\hat{\bm{\sigma}}_{N} \cdot \bm{n}_y  + \gamma \mathrm{d} \bm{u}\right).
$$ (linearized-frictionless-contact)

The case with friction requires some additional consideration.
First, denote $\bm{q} = \hat{\bm{\sigma}}_N - \gamma \dot{\bm{u}}$, $\bm{q}_t = \bm{T}_{n_p} \bm{q}$, and $\delta = -\mathscr{F}\left[ \hat{\sigma}_n(\bm u) + \gamma g(\bm u) \right]_{\mathbb{R}^-}$. Then, the frictional part of the weak form is rewritten as

$$
    -\int_{\Gamma_c} \bm{v}\cdot P_{B\left(\bm{n}_y,\ \delta \right)}(\bm q) \ dS
$$

The linearization of the friction term is then

$$
    \mathrm{d} \left(P_{B(\bm{n}_y, \delta)}(\bm q)\right) =
    \begin{cases}
        0, & \text{if } \delta \leq 0 \text{ or } \|\bm{q}_t\| = 0 \\
        \mathrm{d}\bm{q}_t, &\text{if } \| \bm{q}_t \| \leq \delta \\
        \frac{\delta}{\|\bm{q}_t\|} \left(I - \frac{\bm{q}_t}{\|\bm{q}_t\|} \otimes \frac{\bm{q}_t}{\|\bm{q}_t\|}  \right) \mathrm{d}\bm{q}_t + \frac{\bm{q}_t}{\|\bm{q}_t\|}\mathrm{d}\delta, & \text{otherwise}.
    \end{cases}
$$

Since $\mathrm{d}\delta$ is given in {eq}`linearized-frictionless-contact` and $\mathrm{d} \bm{q}$ can be written as

$$
    \mathrm{d} \bm{q}_t = \mathrm{d} \left( \hat{\bm{\sigma}}_t - \gamma \dot{\bm{u}}_t \right) = \bm{T}_{\bm{n}_y} \left(\mathrm{d} \hat{\bm{\sigma}}_N - \gamma \mathrm{d}\dot{\bm{u}}\right),
$$

the full linearization for contact with friction is given as

$$
    \int_{\Gamma_c} \bm{v}\cdot \mathrm{d}\left(\left[ \hat{\sigma}_n(\bm u) + \gamma g(\bm u) \right]_{\mathbb{R}^-}\right) \bm{n}_y \ dS - \int_{\Gamma_c} \bm{v} \cdot \mathrm{d} \left(P_{B\left(\bm{n}_y, -\mathscr{F}\left[ -\hat{\bm{\sigma}}_N\cdot\bm{n}_y + \gamma g \right]_{\mathbb{R}^-}\right)}(\hat{\bm{\sigma}}_N - \gamma \dot{\bm{u}})\right)\ dS.
$$

## Shapes

Each shape requires an implementation of {c:type}`RatelBCContactProject`, which computes the projection of a given point onto the nearest point $\bm{y}$ on the contact surface and the normal to the contact surface $\bm{n}_y$.
Additionally, all shapes support linearly-interpolated rigid-body translation through space, specified via the `-bc_contact_[face]_translate` option (see {ref}`using-contact` for more details on the command-line interface).
To support these rigid-body motions, each shape requires an implementation of {c:type}`RatelBCContactUpdateModel`.
This function takes a given time and returns an updated shape-specific struct with the geometric properties of the shape at that time.

### Platen
The platen is modeled as an abstract half-space with a point $\bm{c}$ and a normal vector $\bm{n}_p$, which is displaced by a distance $d > 0$ along $\bm{n}_p$ linearly throughout the simulation.
For the platen, finding the nearest point $\bm{y}\in\partial\mathcal{B}$ to a given current configuration point $\bm{x}$ is quite simple.
Note that since the platen has the same normal everywhere, $\bm{n}_y = \bm{n}_p$ and the normal gap can be computed trivially as

$$
    g(\bm{x}) = \bm{n}_p \cdot \left(\bm{x} - \bm{c}\right).
$$

Then, the nearest point is given by subtracting the normal gap from $x$ as:

$$
    \bm{y} = \bm{x} - g(\bm{x}) \bm{n}_p.
$$

For a given time $t$ with prescribed translation $\bm{T}$, the geometric properties are updated as

$$
\begin{align}
\bm{n}_p(t) &= \bm{n}_p\\
\bm{c}(t) &= \bm{c} + \bm{T},
\end{align}
$$

that is, the center $\bm{c}$ is translated.
Platens also support a simplified translation interface via the `-bc_contact_[face]_distance [float list]` option, which prescibes the translation as a distance along the normal vector (see {ref}`using-contact-shape-platen` for more details).