# A novel multi-dimensional *a priori* contact detection algorithm

Sixtine Michel, Gael Poette, Nicolas Therme

► To cite this version:

Sixtine Michel, Gael Poette, Nicolas Therme. A novel multi-dimensional *a priori* contact detection algorithm. 2025. hal-04981784

HAL Id: hal-04981784

https://hal.science/hal-04981784v1

Preprint submitted on 7 Mar 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A Novel Multi-Dimensional *A Priori* Contact Detection Algorithm

Sixtine Michel<sup>a,1</sup>, Gael Poette<sup>b,c,2</sup>, Nicolas Therme<sup>b,3</sup>.

## Abstract

The accurate numerical simulation of moving bodies and interfaces requires robust contact detection algorithms, particularly in high-velocity impact scenarios involving elasto-plastic materials. Traditional contact detection methods, such as node-to-segment (NTS) techniques and neighbor-map approaches, exhibit limitations in both precision and computational efficiency: corrections are typically made after detecting material penetration *a posteriori*, disregarding the dynamics of the configuration. This paper presents a novel multi-dimensional contact detection algorithm based on the Continuous Extension of Thales' Theorem (CETT). The CETT algorithm enables precise determination of contact time and location across one-, two-, and three-dimensional problems. It achieves superior accuracy in *a priori* detection of interface interactions while maintaining computational efficiency, particularly for complex, unstructured meshes. While CETT establishes a new foundation for contact-impact algorithms, it can also be integrated into existing codes. We validate the algorithm through comprehensive numerical benchmarks, demonstrating its robustness and convergence properties in comparison to conventional approaches. The CETT method is fully compatible with both staggered and cell-centered finite volume schemes, making it well-suited for Lagrangian simulations and parallel implementations in industrial and scientific applications.

**Keywords:** Numerical contact detection, Lagrangian framework, euclidean geometry, materials dynamic, Thales theorem.

## 1. Introduction

This paper focuses on the physics of high-velocity impacts between elasto-plastic materials. The Euler equations accurately model the behavior of these materials:

$$ \left\{ \begin{array}{l} \partial_t \rho + \nabla_x \cdot (\rho v) = 0, \\ \partial_t \rho v + \nabla_x \cdot (\rho v \otimes v + \mathbb{T}) = 0, \\ \partial_t \rho e + \nabla_x \cdot (\rho v e + \mathbb{T}^T v) = 0, \\ \mathcal{O}(\mathbb{T}), \end{array} \right. \qquad (1) $$

<footnote><sup>a</sup>University of Ferrara, Via Ariosto, 35, Ferrara, 44121, Italy</footnote>

<footnote><sup>b</sup>CEA CESTA, 15 Avenue des Sablières, Le Barp, 33114, France</footnote>

<footnote><sup>c</sup>Institut de Mathématiques de Bordeaux, Université de Bordeaux, CNRS, Bordeaux INP, France</footnote>

<footnote n="1"><sup>1</sup>Corresponding author. E-mail address: sixtine.michel@unife.it</footnote>

<footnote n="2"><sup>2</sup>Corresponding author. E-mail address: gael-poette@cea.fr</footnote>

<footnote n="3"><sup>3</sup>Corresponding author. E-mail address: nicolas.therme@cea.fr</footnote>
where $x \in \Omega \subset \mathbb{R}^d$ is the spatial variable in the $d$-dimensional spatial domain $\Omega$ and $t \in [0, T] \subset \mathbb{R}^+$ is the time variable. The quantity $\rho = \rho(x, t) \in \mathbb{R}^+$ denotes the mass density, $v = v(x, t) \in \mathbb{R}^d$ represents the material velocity, and $e = e(x, t) \in \mathbb{R}^+$ denotes its total specific energy. System (1) is closed through the Cauchy stress tensor $\mathbb{T}$ and the closure operator $\mathcal{O}$ acting on $\mathbb{T}$.

The closure may depend simply on $\rho$ and $\epsilon = e - \frac{1}{2}\|v\|^2$, yielding $\mathbb{T} = \mathbb{T}(\rho, \epsilon)$, or the operator $\mathcal{O}(\mathbb{T})$ may take the form of an additional partial differential equation [1]. Depending on the structure of operator $\mathcal{O}$ [1], one may consider either a hyperelastic model [2] or a hypoelastic model such as Wilkins' model [3]. Although Wilkins' model is known to be mathematically inconsistent [4], it remains the foundation of many simulation codes. Therefore, in this paper, we develop contact-impact algorithms designed to be minimally sensitive to the specific formulation. Our goal is to handle contact-impact scenarios with Wilkins' model while requiring minimal modifications to existing legacy codes when implementing new, more sophisticated elastoplastic models.

These models, combined with contact-impact simulation capabilities, are highly valuable in industrial and military codes for numerical simulation of high-strain dynamics in solids (see [5, 6, 7, 8] for representative examples). In this context, system (1) is typically solved in Lagrangian coordinates—that is, in the material reference frame. This involves the coordinate transformation:

$$ \left\{ \begin{array}{l} \frac{d}{dt}x(t) = v(x(t), t), \\ x(0) = X, \end{array} \right. \qquad (2) $$

where $v$ represents the material velocity. System (1) in Lagrangian coordinates $X, t$ becomes:

$$ \left\{ \begin{array}{l} \rho_0 \partial_t \tau + \nabla_X \cdot v = 0, \\ \rho_0 \partial_t v + \nabla_X \cdot \mathbb{T} = 0, \\ \rho_0 \partial_t e + \nabla_X \cdot \mathbb{T}^T v = 0, \\ \mathcal{O}_L(\mathbb{T}), \end{array} \right. \qquad (3) $$

where $\tau = \frac{1}{\rho}$ is the specific volume, $\rho_0$ is the initial density, and $\mathcal{O}_L$ is the operator $\mathcal{O}$ expressed in Lagrangian coordinates. While solving (3) is generally more complex than solving (1), it ensures better control of numerical diffusion for contacts (see the seminal works [9, 10] and subsequent studies [11, 12, 13, 14, 15]).

Now, if we come back to our contact-impact problem, the main difficulty in Lagrangian coordinates, comes from the fact that we have as many changes of variables, given by (2), as materials involved in the simulation. And these changes of variables must coincide when the materials are in contact. To get an intuition, it is more convenient to rewrite (3) as an integro-differential equation: for a point $x(t) \in \omega(t)$ (where $\omega(t)$ denotes a moving region of the $d$-dimensional euclidean space) whose displacement is characterized by the trajectory (2), the Lagrangian representation is defined by the motion of the control volume $\omega(t)$ following the fluid velocity. The conservation of the
mass, momentum and total energy read [16]

$$ \left\{ \begin{aligned} & \frac{d}{dt} \int_{\omega(t)} \rho d\nu &&= 0, \\ & \frac{d}{dt} \int_{\omega(t)} \rho v d\nu - \int_{\partial\omega(t)} \mathbb{T}nds &&= 0, \\ & \frac{d}{dt} \int_{\omega(t)} \rho e d\nu - \int_{\partial\omega(t)} \mathbb{T}n \cdot vds &&= 0, \\ & \int_{\omega(t)} \mathcal{O}_L(\mathbb{T}), \end{aligned} \right. \qquad (4) $$

with $\partial\omega(t)$ the boundary of the domain $\omega(t)$, $n$ the outward normal. Now, in Lagrangian coordinates, when several materials are at play, by definition, several reference frames exist and collapse when the materials encounter a contact: take two materials in Lagrangian coordinates, hence two changes of variables similar to (2), one with velocity $v_1$ for material 1 and one with velocity $v_2$ for material 2: then our simulation domain is $\omega(t) = \omega_1(t) \cup \omega_2(t)$ and the boundary is

- $\partial\omega(t) = \partial\omega_1(t) \cup \partial\omega_2(t)$ when there is no contact,

- $\partial\omega(t) = [\partial\omega_1(t) \setminus \partial\omega_1(t) \cap \partial\omega_2(t)] \cup [\partial\omega_2(t) \setminus \partial\omega_1(t) \cap \partial\omega_2(t)] \cup [\partial\omega_1(t) \cap \partial\omega_2(t)]$ when there is a non-empty contact surface $\partial\omega_1(t) \cap \partial\omega_2(t) \neq \emptyset$.

In a sense, the numerical resolution of the contact-impact problem in Lagrangian coordinates can be summed up [17] by

- being able to numerically build $\partial\omega_1(t) \cap \partial\omega_2(t)$ when it is different from $\emptyset$. This step is called *detection*. It is knownly complex [18, 19]. Some strategies consequently aim at astutely avoiding it, see [20, 21]. For these strategies, this step and the next one are merged.

- Being able to define the changes of variables for the different materials in $\omega(t) = \omega_1(t) \cup \omega_2(t)$ together with having them to coincide in $\partial\omega_1(t) \cap \partial\omega_2(t) \neq \emptyset$. This step is commonly called the *correction*.

Several strategies are suggested in the literature, for both of these steps or in order to make them quite transparent. These will be discussed later on. But before, we suggest detailing a method, implemented in a 1D legacy code of the CEA$^4$, which to our knowledge, has never been published and which presents interesting properties. It will represent the starting point of the study documented in this paper: let us consider a simple 1D Taylor impact in a Lagrangian context. This type of experiments is intensively used in order to characterize parameters involved in some material models for $\mathbb{T}$ for example. All the values are dimensionless for the sake of simplicity. Assume we initially have two materials $\omega_1(0) = [x_1^l(0), x_1^r(0)] = [1, 2]$ and $\omega_2(t) = [x_2^l(0), x_2^r(0)] = [2, 4]$. Both materials are initially in contact and are the same (see Fig. 1): their closure is modeled by a Mie-Gruneisen equation of state, i.e. $\mathbb{T} = p$ with reference density $\rho_0^r = 10$, bulk

<footnote>$^4$The authors of this paper do not claim the paternity of 1D contact method described in this introduction. The idea comes from J. M Favier (retired from the CEA)</footnote>
<figure><img src="image_2.png" /><figcaption>Figure 1: Initial condition for the one-dimensional contact problem of the introduction.</figcaption></figure>

modulus $\Gamma = 2$ and Gruneisen coefficient $K_s = 10$. Initially, the second material is at rest ($v_2 = 0$) and the first one has velocity $v_1 = (v_x^1, v_y^1, v_z^1)^T = (1, 0, 0)^T$. This means the material are going to impact at time $t^* = 0$. Otherwise, material 1 and 2 are the same in the sense that $\rho_1(0) = 10 = \rho_2(0)$ and $p_1(0) = p_2(0) = 0$. As soon as $t > 0$, an impact occurs and some waves propagate through the materials. These are such that after a certain time both materials separate. But, due to the left hand side boundary condition which is an imposed pressure of $p(x_1^l(t), t) = 5$, material 1 comes back toward material 2 and they are back in contact until final time $T = 8$ is reached. Figure 2 (left) gives an idea of the dynamics: the curves $t \mapsto \omega_1(t)$ and $t \mapsto \omega_2(t)$ are displayed for several meshes with $N_x \in \{150, 300, 600, 1200, 2400, 4800, 9600, 19200, 38400\}$ cells. The material remains in contact for $t \lesssim 3.5$. At $t \approx 3.5$, material 2 separates from material 1. But as the left hand side boundary condition on material 1 is an imposed pressure law, material 1 is pushed against material 2 and contact is made at $t \approx 5.9$. The dots on figure 2 (left) correspond to the impact position and time $(x^*, t^*)$, i.e. $t^* = \inf\{t > 0 | x_1^r(t) = x_2^l(t)\}$ and $x_1^r(t^*) = x_2^l(t^*) = x^*$, for each spatial discretization $N_x \in \{150, 300, 600, 1200, 2400, 4800, 9600, 19200, 38400\}$. This simple test-case allows

<figure><img src="image_3.png" /><figcaption>Figure 2: Left: curves $t \mapsto \omega_1(t)$ and $t \mapsto \omega_2(t)$ of the boundaries of materials 1 and 2 which may be in contact for several meshes. Right: curves $t \mapsto v_2(x_2^r(t), t)$ for the free surface velocity with respect to time for several meshes.</figcaption></figure>

exploring several crucial phases of the contact-impact physics: material separation, detection and correction.

Now, in an industrial context, when dealing with such kind of set-up, we are usually interested in being able to recover the free surface velocity of the second material with respect to time, i.e. $t \mapsto v_2(x_2^r(t), t)$. This quantity can be measured in experiments [22, 23, 24, 25, 26, 27], the simulation code and the models within can be validated. Figure 2 (right) displays the free surface velocity profiles for meshes with $N_x \in \{150, 300, 600, 1200, 2400, 4800, 9600, 19200, 38400\}$ cells: initially, the velocity is
zero. But the shock, initiated by the difference of initial velocities in the material, goes through material 2 and makes the free right hand side surface of material 2 move: its velocity jumps to 2.0 and oscillates between 2.0, 2.6 and 0.5 until time $t^* \approx 5.9$ which corresponds to the time at which material 1 reaches material 2. After this reshock, the free surface velocity is once again accelerated.

On both pictures of figure 2, qualitative converging behaviors are observable testifying on the accuracy of both the numerical scheme implemented in the code (Godunov scheme [28] here) and the detection-correction strategy. Figure 3 is more quantitative

<figure><img src="image_4.png" /><figcaption>Figure 3: Top: convergence study of the impact position and time ($x^*, t^*$). An order 1 with respect to $\Delta x$ converging behaviour is observed for both observables. Bottom: convergence study for the free surface velocity $v_2(x_2^r(t = 7.5), t = 7.5)$. Same observation: an order 1 convergence is observed as soon as a certain mesh size is reached.</figcaption></figure>

with respect to convergence: it displays some convergence studies for three physical observables, namely the impact time $t^* = \inf\{t > 0|x_1^r(t) = x_2^l(t)\}$, the impact position $x^* = x_1^r(t^*) = x_2^l(t^*)$, the free surface velocity $t \mapsto v_2(x_2^r(t = 7.5), t = 7.5)$. From figure 3, we observe an $O(\Delta x)$ convergence for every observables. This property is mandatory in order to avoid error compensations during calibration for example. This ensures being able to choose a mesh such that the discretisation error will be negligible with respect to every other sources of errors/uncertainties. In practice, we obtain the same curves with the VNR [29] scheme (i.e a staggered scheme instead of a centered one), with the same contact-impact strategy.
From these converging behaviours, we infer that an $\mathcal{O}(\Delta x)$ contact-impact strategy will ensure an $\mathcal{O}(\Delta x)$ convergence on the post impact observables of interest.

Now, in practice, on the previous benchmark, in order to obtain these $\mathcal{O}(\Delta x)$ converging behaviours for the contact-impact strategy for the centered and staggered schemes,

- detection is a priori handled by numerically looking for time $t^* = \inf\{t > 0 | x_1^r(t) = x_2^l(t)\}$. In practice, this is done thanks to the fact that on a given time step $[t^n, t^{n+1} + \Delta t = t^{n+1}]$, we have access to the velocities $v_1^n, v_2^n$ of each boundary node of the mesh. We consequently have $\forall t \in [t^n, t^{n+1}], \alpha \in \{1, 2\}$, $x_\alpha(t) = u_\alpha^n(t - t^n) + x_\alpha^n$. And $t^*$, defined by $x_1^r(t^*) = x_2^l(t^*)$ is the solution, if it exists, of a first order equation. In practice, if $t^* \in [t^n, t^{n+1} = t^n + \Delta t]$, the time step is decreased to $\min(\Delta t, t^* - t^n)$ in order to reach $t^*$ at the end of the iteration (if contact occurs).

- The correction now is only based on the fact that the different nodes are dealt as boundary conditions or not. Basically, the connectivity table of the cells changes depending on whether the flux between two cells must be computed from a neighboring cell (contact) or from a boundary condition (separation). This phase may slightly differ depending on whether a centered or a staggered scheme is considered:

- for a centered scheme, the flux between the cells in contact is classically computed. One only needs to flag the nodes, they are either free or in contact and the two points must be considered respectively at the boundary or within the domain.

- For a staggered one, the nodes $x_1^r(t^*)$ and $x_2^l(t^*)$ must be such that the change of variable between the two materials is unique at this point. In practice, conservation of the nodal mass, nodal momentum and nodal total energy is enforced.

From the description of the two phases, first, note that for the detection, we never really needed to

- neither define the numerical scheme used in the code (we just need the node velocities in order to extract $t^*$),

- nor the model we solve (hyperelasticity, Wilkins), we rely on continuous principles (conservation).

For the correction, the discussions, slightly different for a centered scheme and a staggered one, are closely related to the question *how to deal with non conform meshes?*.

In other words, the contact-impact algorithm previously described only need boundary nodal velocities (even if accelerations are used if available). If those quantities can accurately enough be estimated from the numerical scheme solving a given model, then the previous detection-correction steps remain independent of the solver and of the model. This property is important: this means that the implemented contact algorithm remains the same whether a staggered or a centered scheme is used (and
both can be used with contact if they are implemented in the same code), to solve an
hyperelastic system or an hypoelastic one. The correction step may be slightly differ-
ent depending on the kind of scheme (just as boundary conditions are usually treated
differently) but the general architecture of the algorithm remains the same. To our
knowledge, these properties are quite singular in the literature as there are two main
families of contact treatments:

- *Variational algorithms:* these methods originated in structural mechanics ap-
plications requiring high precision in contact conditions, including a strict non-
penetration constraint. Initially developed for implicit finite element codes, La-
grangian contact methods were designed to solve large deformation and quasi-
static problems by enforcing impenetrability constraints practically ensured by
solving an optimization problem [30, 31, 32, 33, 34]. These methods have been
the subject of more recent developments [21, 35, 36] including an adaptation
within the context of centered finite volume scheme for Lagrangian hydrodynam-
ics [20, 37]. The main idea is to position the nodes of the material by solving
a minimisation problem of size $N_x \times N_x$ with impenetrability constraints ($N_x$ is
the total number of cells in every materials). Note that the term *Lagrangian*
in the previous denomination does not refer to the Lagrangian coordinates but
rather to the Lagrangian multiplicators used in order to treat the constraints of
the optimization problem. The constraints are quite tricky to define [20], the
minimisation problem may be costly but convergence is guaranteed [20, 37].

- *A posteriori contact detection-correction algorithms:* These algorithms have been developed to handle contact in fast transient dynamics codes [8, 5, 7]. These codes are often explicit in time, and the methods mentioned previously are difficult to apply because of their cost. In this type of algorithm, each material evolves freely at each iteration, and interpenetrations are detected and corrected *a posteriori*. The main idea is to provide a local treatment on a reduced stencil to avoid coupling the degrees of freedom of the problem. The correction part consists in modifying the interpenetrating nodes position velocity and acceleration in order to impose contact conditions. The methods commonly used are the penalty method [38], explicit Lagrange multipliers (displacement constraint imposed at the next iteration) [39, 40], and Lagrange multipliers [41]. The interested reader can refer to [42] for more in depth description of such methods. These corrections are generally cheap but their local nature, and the approximations introduced, tend to create mesh sensitivity which can make convergence more difficult (user tuning is sometimes necessary).

As for the detection procedure, the classical Node To Segment approach (NTS)
[18] is widely used in explicit legacy hydrodynamics code and for a long time now
[42]. This algorithm consists in determining the penetration of the nodes of a
"slave" surface through the faces of a "master" surface. Variations exist on this
algorithm, mainly by simplifying the contact geometry to speed up calculations,
for example by using spheres to represent the contact surfaces [43, 44]. For more
details, a large overview of contact detection techniques in hydrocodes can be
found in [17]. We also note that contact detection techniques are widely devel-
oped and used in video games for decades now (see [45] and references therein),
however they are not used for hydrodynamic problems which need high accuracy.
The NTS algorithm will be considered the baseline for detection comparison
throughout this paper.

This paper begins with the observation that, to our knowledge, none of the previous
2D/3D contact strategies described in the literature can be understood as a general-
isation of the (satisfactory) 1D one described in this introduction. We consequently
quite naturally try to extend it, in this paper, to 2D and 3D problems. The aim of
this paper is to document how the detection phase can be efficiently handle in 2D and
3D, with an algorithm which naturally degenerates toward the one described in this
introduction for 1D problems, ensuring an $O(\Delta x)$ convergence. We will see that the
detection of contact, quite simple in 1D, deserves some special attention in 2D/3D.
Although there is a large literature on the subject, there is little work comparing de-
tection algorithms. In this paper we will propose a comparison between two of them:
the first one is the previously mentionned NTS technique (detailed in §4.1), considered
as a baseline for detection, and the second one is the purpose of this paper, is based
on a kind of generalisation of Thalès' theorem. In this paper, we aim at designing a
robust and accurate contact detection procedure fully parallelizable, and usable with
staggered, but also centered schemes, independent of the (hyper-hypo) elasto-plastic
model.

The paper is organised as follows: in §2, we introduce notations and explicit the
contact detection procedure from a continuous point of view. In §3, we discretize
the continuous problem, show that the new contact detection algorithm we suggest
naturally degenerates to the one previously described in 1D and explicit the differences
and the numerical choices for the 2D/3D cases. We describe the NTS procedure (our
baseline for comparisons) in §4.1. In particular, we highlight exceptions and singular
cases where this algorithm is not reliable or deserves particular attention. Finally, §4 is
devoted to numerical benchmarks. In particular, care is taken to verify our algorithm
exhibits the expected convergence rates.

## 2. The Contact Detection Procedure: A Continuous Formulation

This section presents a continuous formalization of *contact detection*. The multidimensional discretizations presented in §3 (1D/2D/3D) are all derived from the framework established here.

Consider two moving bodies in an *n*-dimensional space, as illustrated in Figure 4.
These bodies are defined by $x \in \omega_1$ and $x \in \omega_2$, with their respective interfaces denoted
by $\partial\omega_1$ and $\partial\omega_2$.
<figure><img src="image_5.png" /><figcaption>Figure 4: Two closed spaces $\omega_1$ and $\omega_2$.</figcaption></figure>

We assume that the displacement of both bodies is fully defined by the mappings $t \mapsto \omega_1(t)$ and $t \mapsto \omega_2(t)$. The detection phase aims to determine whether contact or intersection occurs during the time interval $t \in [0, T]$ and, if so, to precisely identify both the time of occurrence within $[0, T]$ and the location within the physical space $\omega_1(t) \cup \omega_2(t)$.

Computational contact detection has been extensively studied in the literature. As discussed in Section 1, the treatment of contact through specialized algorithms is crucial in Lagrangian coordinates: the materials are naturally separated entities, each defined solely through their respective coordinate transformations of the form (2). This separation means that Lagrangian materials do not inherently interact with each other. Consequently, in Lagrangian simulations involving multiple materials, inter-penetrations may occur without proper detection, leading to non-physical computations (see Figure 5).

<figure><img src="image_6.png" /><figcaption>Figure 5: Illustration of potential inter-penetration: A static body $\omega_1$ and a moving body $\omega_2$. Without a contact-detection-correction algorithm, $\omega_2$ moves through $\omega_1$ during the simulation, ignoring its presence.</figcaption></figure>

To define contact between the two entities $t \mapsto \omega_1(t)$ and $t \mapsto \omega_2(t)$, we require knowledge of their mappings, which depends on the velocities of all material points within these entities. Specifically, for any position $t \mapsto x_1(t)$ and $t \mapsto x_2(t)$ of materials 1 and 2, we must have access to their respective velocity mappings $t \mapsto v_1(t) = v_1(x_1(t))$ and $t \mapsto v_2(t) = v_2(x_2(t))$, as shown in Figure 6.
<figure><img src="image_7.png" /><figcaption>Figure 6: Representation of two nodes $x_1(t) \in \partial\omega_1(t)$ and $x_2(t) \in \partial\omega_2(t)$, and their respective velocity fields $v_1(x_1(t))$ and $v_2(x_2(t))$.</figcaption></figure>

Let us consider we are in an Euclidean space of dimension $n \in \mathbb{N}$. The norm attached to this space is the $L_2$-norm denoted by $||.||$. Let us now define a contact: a contact will occur at the time $t^* \in [0, T] \subset \mathbb{R}^+$ if and only if $\exists t^* \in [0, T]$ such that $\partial\omega_1(t^*) \cap \partial\omega_2(t^*) \neq \emptyset$.

In time interval $[0, t]$, the trajectory of any material point $x(t)$, independently of being from one material or another, is given by

$$t \mapsto x(t) = x(0) + \int_{0}^{t} v(s)ds = x_{0} + \int_{0}^{t} v(x(s))ds.$$

Let us now consider $x_1(t) \in \partial\omega_1(t)$ and $x_2(t) \in \partial\omega_2(t)$, two points on the boundaries of materials 1 and 2. Let us define $t^*$ by

$$t^* = \inf\{t > 0 | \exists (x_1(t), x_2(t)) \in \partial\omega_1(t) \times \partial\omega_2(t) | x_1(t) = x_2(t)\}.$$

This time of contact is defined by the existence of two material points $x_1, x_2$ of the surfaces $\partial\omega_1, \partial\omega_2$ coinciding for the first time at $t^*$. With the previous definition, $t^*$ is quadratic in terms of complexity as every position on the surface of material 1 must be compared to every position on the surface of material 2. The matter of §3 will be to alleviate the more possible this complexity. Of course, the problem will remain quadratic. But we will make sure it is tractable. Contact occurs in interval of time $[0, t]$ if $t^* \in [0, t]$. The same $t^*$ can be expressed with respect to the velocities and initial conditions as we have

$$t^* = \inf \left\{ t > 0 \middle| \exists (x_1(0), x_2(0)) \in \partial\omega_1(0) \times \partial\omega_2(0) | x_1(0) + \int_0^t v_1(s)ds = x_2(0) + \int_0^t v_2(s)ds \right\}.$$

This means that in practice, it is enough having the positions on the surfaces of the materials and their respective velocity fields. Let us work even more on $t^*$ and rewrite the above expressions as

$$t^* = \inf \left\{ t > 0 \middle| \exists (x_1(0), x_2(0)) \in \partial\omega_1(0) \times \partial\omega_2(0) | x_1(0) - x_2(0) = \int_0^t [v_2(s) - v_1(s)] ds \right\}. \quad (5)$$

With the previous expression, the existence of $t^*$ (in the sense $0 < t^* < \infty$) is closely related to the colinearity of $x_1(0) - x_2(0)$ and $v_2(s) - v_1(s)$, $\forall s \in \mathbb{R}^+$ during the trajectory. In practice, in the following sections, we will calculate $t^* \in \mathbb{R}^+$ and a posteriori
check whether it is within $[0, t]$ to detect contact or not. The contact locations are quite intuitively defined via $t^*$ as the points $x_1(t^*) = x_2(t^*) = x(t^*) \in \partial\omega_1(t^*) \cap \omega_2(t^*)$.

In the next section, from the previous continuous definition, we suggest discretisations which are "solver aware" in the sense that these discretisations must be compatible with the discretisation schemes of system (3) (independently of being centered of staggered) and suggest ways to alleviate the complexity needed to solve a quadratic problem with respect to all the points of the boundaries.

### 3. Multi-dimensional Discretization of the CETT

In this section, we extend our previous continuous framework to a discrete formulation that aligns with numerical solver requirements. Specifically, we examine the implementation across dimensions from 1D to 3D. In Section §3.1, we demonstrate how our detection strategy naturally reduces to the approach described in §1. The 2D/3D sections then present a generalization of this 1D procedure, with details developed in Sections §3.2–§3.3. The new procedure is called CETT (Continuous Extension of Thales’ Theorem), a name that will be justified in Section §3.2. Throughout these sections, we focus on developing computationally efficient methods to find $t^*$ while maintaining $O(\Delta x)$ convergence order.

#### 3.1. The One-dimensional Contact Detection Algorithm

In one-dimensional configurations, the discretization of (5) is straightforward since boundaries in 1D are simply nodes (see Fig. 7).

<figure><img src="image_8.png" /><figcaption>Figure 7: One-dimensional problem. Two materials are represented. Their respective boundaries are the extremal nodes.</figcaption></figure>

Let us consider material 1 on the left-hand side and material 2 on the right-hand side. We denote the boundary nodes of materials 1 and 2 by $\partial\omega_1(t) = \{x_1^r(t), x_1^l(t)\}$ and $\partial\omega_2(t) = \{x_2^r(t), x_2^l(t)\}$ respectively. In 1D, each material’s boundary is fully defined by two nodes. Given the a priori knowledge that material 1 is on the left-hand side of material 2, the continuous definition (5) becomes

$$t^* = \inf \left\{ t > 0 | x_1^r(0) - x_2^l(0) = \int_0^t [v_2^l(s) - v_1^r(s)] ds \right\}. \quad (6)$$

The quadratic problem reduces to a single equation. At this stage, equation (6) remains continuous. Assuming the time interval $[0, t]$ corresponds to a time step of a resolution scheme solving (3), the velocities remain constant during $[0, t]$. Equation (6) then becomes

$$ \begin{aligned} t^* &= \inf \{t > 0 | x_1^r(0) - x_2^l(0) = t [v_2^l(0) - v_1^r(0)] \}, \\ &= \frac{v_2^l(0) - v_1^r(0)}{x_1^r(0) - x_2^l(0)}. \end{aligned} \quad (7) $$
The existence of $t^*$ is closely related to having collinear vectors of initial conditions for positions and velocities. Most importantly, the discretization of (5) naturally leads to the numerical strategy described in §1. Thus, the 2D and 3D descriptions that follow are generalizations of what we have just presented.

**Remark:** The numerical method described in §1 is not exactly the discretization of expression (5). More precisely, the discretization of (5) is

$$t^* = \min \left\{ \frac{v_2^l(0) - v_1^r(0)}{x_1^r(0) - x_2^l(0)}, \frac{v_2^l(0) - v_1^l(0)}{x_1^l(0) - x_2^l(0)}, \frac{v_2^r(0) - v_1^r(0)}{x_1^r(0) - x_2^r(0)}, \frac{v_2^r(0) - v_1^l(0)}{x_1^l(0) - x_2^r(0)} \right\},$$

where all possible boundary combinations are tested and the minimum is selected. However, some a priori knowledge (specifically, that only certain boundary points can interact) is used to reduce the problem's complexity. Such a priori knowledge is valuable, and our proposed strategy can naturally incorporate it.

Consider $M$ materials that may impact, each with boundaries defined by two nodes. The complexity of calculating $t^*$ at each time step is $O(2^M)$ with an additional min computation. When a priori knowledge is incorporated (i.e., materials are ordered), this complexity reduces to $O(M - 1)$ plus a min operation.

Moving forward, we assume that materials are discretized using polygonal meshes, which is typical when using finite volumes or low-order finite elements. The method we present in the following sections for 2D and 3D assumes straight edges and faces that can be decomposed into planar sub-faces.

## 3.2. The two-dimensional contact detection algorithm

In this section, we follow the same methodological approach as in the previous section: first, we define the discretized surfaces $\partial\omega_1, \partial\omega_2$, then revisit the expression (5) of $t^*$ in this context and introduce the discretization scheme hypotheses. Finally, we demonstrate that the problem reduces to detecting contact between an arbitrary node and an arbitrary face of the two surfaces, and we present an efficient computational method.

For the two-dimensional problem, to define the discretization of $t^*$, we must first define the boundary surfaces of $\omega_1(t)$ and $\omega_2(t)$. The material surfaces are defined through their respective nodes $(x_{\beta,i})_{\beta \in \{1,2\},i \in \{1,...,N_{\beta}\}}$, where $N_{\beta}$ is the number of nodes in $\partial\omega_{\beta}$, the surface of material $\beta \in \{1,2\}$. In 2D, we can assume that the nodes are ordered such that $f_{i,i+1}^{\beta} = [x_{\beta,i}, x_{\beta,i+1}]$ defines a face for material $\beta$ between two adjacent nodes $x_{\beta,i}$ and $x_{\beta,i+1}$. This is illustrated in figure 8. Consequently, the material surfaces can be defined equivalently either through their nodes $\partial\omega_{\beta} = \cup_{i \in \{1,...,N_{\beta}\}} x_{\beta,i} = \{x_{\beta,1}, ..., x_{\beta,N_{\beta}}\}$ or through their faces $\partial\omega_{\beta} = \cup_{i \in \{1,...,N_{\beta}\}} f_{i,i+1}^{\beta} = f_{1,2}^{\beta} \cup ... \cup f_{N_{\beta}-1,N_{\beta}}^{\beta} = [x_{\beta,1}, x_{\beta,2}] \cup ... \cup [x_{\beta,N_{\beta}-1}, x_{\beta,N_{\beta}}]$. This dual representation of discrete surfaces will be extensively used in the following algorithms.

Now, let us revisit expression (5) with our 2D discretization. The first challenge in 2D arises from the fact that $\partial\omega_1$ and $\partial\omega_2$ are not discrete: in expression (5), we must consider all initial points on the boundaries of both materials. To construct a finite-dimensional problem, where the min is taken over a finite number of computed times, we leverage the dual representation of the surfaces. Indeed, expression (5), recalled
<figure><img src="image_9.png" /><figcaption>Figure 8: Two dimensional problem: materials 1 and 2 are defined via their nodes $(x_{1,j})_{j \in \{1..6\}}$ for $\partial\omega_1$ and $(x_{2,j})_{j \in \{1..6\}}$ for $\partial\omega_2$. For material $\beta \in \{1,2\}$, the faces between two nodes $x_{\beta,i}$ and $x_{\beta,i+1}$ is denoted by $f_{i,i+1}^{\beta} = [x_{\beta,i}, x_{\beta,i+1}]$.</figcaption></figure>

below

$$t^* = \inf \left\{ t > 0 \middle| \exists (x_1(0), x_2(0)) \in \partial\omega_1(0) \times \partial\omega_2(0) \middle| x_1(0) - x_2(0) = \int_0^t [v_2(s) - v_1(s)] ds \right\}, \quad (8)$$

is equivalent in 2D to

$$t^* = \inf \left\{ \begin{array}{l} t > 0 | \exists x_1(0) \in \cup_{i=1}^{N_1} x_{1,i}, x_2(0) \in \cup_{i=1}^{N_2} f_{i,i+1}^2 \\ \text{such that } x_1(0) - x_2(0) = \int_0^t [v_2(s) - v_1(s)] ds \end{array} \right\}, \quad (9)$$

where we have chosen to iterate through a list of points for material 1 and a list of faces for material 2. Note that the reverse approach would be equally valid.

**Remark:** In the literature, surfaces are often categorized as either 'master' or 'slave'. In our context, the master surface is defined by its faces, while the slave surface is defined by its points. The algorithm is not symmetric - swapping the roles of master and slave may yield different results. For this reason, most contact-impact algorithms of this type typically perform two passes, where each material alternates between master and slave roles. This approach symmetrizes the treatment of contact nodes.

Using this expression for $t^*$, the CETT procedure must evaluate each possible combination of a node $x_1$ from material 1 and a face $f_{i,i+1}^2$ from material 2. Let us consider one such point-face pair: at this stage, each face still represents an infinite set of points. This configuration is illustrated in Figure 9. To reduce this to a discrete case, we detect penetration of $x_1$ through a reparametrization of $f_{i,i+1}^2$. A point $x_\alpha$ belongs to $f_{i,i+1}^2$ if there exists $\alpha \in [0, 1]$ such that $x_\alpha = \alpha x_{2,i} + (1-\alpha)x_{2,i+1}$, as shown in Figure 10. Note that this system has its own velocity fields $v_1$, $v_{2,i}$ and $v_{2,i+1}$ (respectively). Using this linear combination, we can define $v_\alpha = \alpha v_{2,i} + (1-\alpha)v_{2,i+1}$.
<figure><img src="image_10.png" /><figcaption>Figure 9: Example of one combination of node $x_1 \in \partial\omega_1$ and surface $f_{i,i+1}^2 = [x_{2,i}, x_{2,i+1}] \in \partial\omega_2$ for which contact detection is tested.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 10: Example of the node $x_1 \in \partial\omega_1$ and the surface $f_{i,i+1}^2 = [x_{2,i}, x_{2,i+1}] \in \partial\omega_2$. Contact detection with the node $x_\alpha \in f_{i,i+1}^2 \in \partial\omega_2$.</figcaption></figure>

Now, suppose that during interval of time $[0, t^*]$, the velocities are constant, then the colinearity condition between $x_1(0) - x_\alpha(0)$ and $v_1(s) - v_\alpha(s), \forall s \in [0, t^*]$ given by

$$x_1(0) - x_\alpha(0) = \int_0^{t^*} [v_\alpha(s) - v_1(s)] ds$$

becomes

$$x_1(0) - x_\alpha(0) = t^* [v_\alpha(0) - v_1(0)]. \qquad (10)$$

Colinearity means we are looking for an $\alpha \in [0, 1]$ such that

$$\overrightarrow{x_1(0)x_\alpha(0)} \wedge (\overrightarrow{x_1(0)} + t^*v_1(0))(\overrightarrow{x_\alpha(0)} + t^*v_\alpha(0)) = 0.$$

This problem is equivalent to finding an $\alpha \in [0, 1]$ such that $\Gamma(\alpha) = 0$, where $\Gamma \in \mathbb{P}_2$ (the set of second order polynomials) defined by$^5$

$$\begin{aligned} \Gamma_{1,2}(\alpha) &= \frac{1}{t^*} \overrightarrow{x_1 x_\alpha} \wedge (\overrightarrow{x_1 + t^* v_1})(\overrightarrow{x_\alpha + t^* v_\alpha}), \\ &= \frac{1}{t^*} (\overrightarrow{x_1 x_\alpha} \wedge \overrightarrow{x_1 x_\alpha} - t^* \overrightarrow{x_1 x_\alpha} \wedge \overrightarrow{v_\alpha v_1}), \\ &= (x_1 - x_\alpha) \wedge (v_\alpha - v_1) = (0, 0, \gamma_z(\alpha)). \end{aligned} \qquad (11)$$

<footnote>$^5$we drop the time dependence (0) for the sake of conciseness.</footnote>
Then, using the definition of $x_\alpha = \alpha x_{2,i} + (1 - \alpha)x_{2,i+1}$ in the above expression leads to

$$ \Gamma_{1,2}(\alpha) = \alpha^2 (\overrightarrow{x_{2,i}x_{2,i+1}} \wedge \overrightarrow{(v_{2,i+1} - v_{2,i})}) \\ + \alpha (\overrightarrow{x_1x_{2,i+1}} \wedge \overrightarrow{(v_{2,i} - v_{2,i+1})} + \overrightarrow{x_{2,i+1}x_{2,i}} \wedge \overrightarrow{(v_{2,i+1} - v_1)}) \\ + (\overrightarrow{x_1x_{2,i+1}} \wedge \overrightarrow{(v_{2,i+1} - v_1)}), \qquad (12) \\ = a\alpha^2 + b\alpha + c. $$

The coefficients $a, b, c$ are introduced mainly for the sake of conciseness and are used in Algorithm 1 later on.

Let us denote by $\alpha^* \in [0, 1]$ the value such that $\Gamma_{1,2}(\alpha^*) = 0$. Assuming $\alpha^*$ exists and is unique, we can compute a contact time $t_c^{1,2}$ between the two points (hence between point $x_1$ and face $f^2$). We can then determine $t_c^{1,2}$ and the contact point $x_\alpha(t_c^{1,2}) = x_1(t_c^{1,2})$. If $t_c^{1,2} < \Delta t$, contact will occur during the next time step between point $x_1$ and face $f^2$. Note that $t_c^{1,2}$ may not equal $t^*$: there may exist another point $x_i/\text{face } f^j$ pair with a smaller contact time $t_c^{i,j} < t_c^{1,2}$. The contact time $t^*$ is the minimum among all possibilities.

In what follows, we detail the computation of $t_c^{1,2}$. For notational simplicity, we omit the superscript $1,2$ in the following paragraph.

We now compute $t_c$. The contact time $t_c$ between $x_1$ and $f_2$ is obtained by computing the scalar product between $x_1(0) - x_{\alpha^*}(0)$ and

$$ x_1(0) - x_{\alpha^*}(0) = t^*[v_{\alpha^*}(0) - v_1(0)]. \qquad (13) $$

This leads to

$$ t^* = \frac{\langle x_1(0) - x_{\alpha^*}(0), x_1(0) - x_{\alpha^*}(0) \rangle}{\langle v_{\alpha^*}(0) - v_1(0), x_1(0) - x_{\alpha^*}(0) \rangle}. \qquad (14) $$

**Remark:** Now, due to the colinearity of vectors $\overrightarrow{x_{1,j}x_\alpha}$ and $\overrightarrow{(x_{1,j} + v_{1,j})(x_\alpha + v_\alpha)}$ in Fig. 10, we know (from Thalès' theorem) that $\frac{\|\overrightarrow{Ox_{1,j}}\|}{\|O(x_{1,j} + v_{1,j})\|} = \frac{\|\overrightarrow{Ox_\alpha}\|}{\|O(x_\alpha + v_\alpha)\|} = \frac{\|\overrightarrow{x_{1,j}x_\alpha}\|}{\|(x_{1,j} + v_{1,j})(x_\alpha + v_\alpha)\|}$. In a sense, the detection procedure we just described, in 2D, is mainly an application of Thalès' theorem. Note also that it is quite easy checking that the definition (14) of $t^*$ is in agreement with the 1D one of (7). This is why the strategy we suggest in this paper is called CETT for continuous extension of Thalès' theorem.

To go further in the resolution of $\Gamma(\alpha) = 0$, two cases can be highlighted:

* first, let us consider the case where $a, b, c$ are not all zero. Depending on the discriminant sign of (12), zero, one or two solutions of $\Gamma_{1,2}(\alpha) = 0$ exist. If $\alpha$ is solution of the equation and $\alpha < 0$ or $\alpha > 1$, it means the point of contact is out of the segment $f_{i,i+1}^2$ and therefore contact is not detected. Finally, root(s)
of $\Gamma_{1,2}(\alpha)$ in $[0, 1]$ are candidates to define the first point of contact between the two domains. The selected one is the one for which the contact time $t_c$ is the smallest.

* Then, in the other case, if the three coefficients of $\Gamma_{1,2}$ (i.e. $a, b, c$) are zero, then $\forall \alpha \in [0, 1]$, $\exists t_{\alpha}^* \in \mathbb{R}$ such as $x_1(t) + (t_{\alpha}^* - t)v_1(t) = x_{\alpha}(t) + (t_{\alpha}^* - t)v_{\alpha}(t)$. By definition of the coefficients $a, b, c$, we have

$$ 
\left\{
\begin{array}{l}
a = 0 \Leftrightarrow \overline{x_{2,i}x_{2,i+1}} \wedge \overline{(v_{2,i+1} - v_{2,i})} = 0, \\
c = 0 \Leftrightarrow \overline{x_1x_{2,i+1}} \wedge \overline{(v_{2,i+1} - v_1)} = 0, \\
b = 0 \Leftrightarrow \underbrace{\overline{x_1x_{2,i+1}} \wedge \overline{(v_{2,i} - v_{2,i+1})}}_{\substack{\text{from } c=0 \\ \overline{x_1x_{2,i+1}} \wedge \overline{(v_{2,i}-v_1)}}} + \underbrace{\overline{x_{2,i+1}x_{2,i}} \wedge \overline{(v_{2,i+1} - v_1)}}_{\substack{\text{from } a=0 \\ \overline{x_{2,i+1}x_{2,i}} \wedge \overline{(v_{2,i}-v_1)}}} = 0, \\
\Rightarrow (\overline{x_1x_{2,i+1}} + \overline{x_{2,i+1}x_{2,i}}) \wedge \overline{(v_{2,i} - v_1)} = \overline{x_1x_{2,i}} \wedge \overline{(v_{2,i} - v_1)} = 0.
\end{array}
\right.
\quad (15) $$

The direct consequence of having $c = b = 0$ is that $\exists \mathcal{I} \in \mathbb{R}, \mathcal{I}' \in \mathbb{R}$ such that (colinearity conditions)

$$ 
\left\{
\begin{array}{l}
\mathcal{I}' \overline{x_1 x_{2,i+1}} = \overline{(v_{2,i+1} - v_1)}, \quad \mathcal{I}' = \frac{\|\overline{(v_{2,i+1} - v_1)}\|}{\|\overline{x_1 x_{2,i+1}}\|}, \quad \left\langle \frac{\overline{x_1 x_{2,i+1}}}{\|\overline{x_1 x_{2,i+1}}\|}, \frac{\overline{(v_{2,i+1} - v_1)}}{\|\overline{(v_{2,i+1} - v_1)}\|} \right\rangle, \\
\mathcal{I} \overline{x_1 x_{2,i}} = \overline{(v_{2,i} - v_1)}, \quad \mathcal{I} = \frac{\|\overline{(v_{2,i} - v_1)}\|}{\|\overline{x_1 x_{2,i}}\|} \quad \left\langle \frac{\overline{x_1 x_{2,i}}}{\|\overline{x_1 x_{2,i}}\|}, \frac{\overline{(v_{2,i} - v_1)}}{\|\overline{(v_{2,i} - v_1)}\|} \right\rangle.
\end{array}
\right.
\quad (16) $$

NB: Rewriting the time of contact $t^*$ (14) as follow

$$ t^*(\alpha) = \frac{\|\overline{x_1 x_{\alpha}}\|}{\|(v_{\alpha} - v_1)\|} \left\langle \frac{\overline{x_1 x_{\alpha}}}{\|\overline{x_1 x_{\alpha}}\|}, \frac{\overline{(v_{\alpha} - v_1)}}{\|(v_{\alpha} - v_1)\|} \right\rangle, \quad (17) $$

if $v_{2,i} \neq v_1$, $\mathcal{I} = t^*(\alpha = 1)$ and if $v_{2,i+1} \neq v_1$, $\mathcal{I}' = t^*(\alpha = 0)$. In this sense, $\mathcal{I} \le 0$ (respectively $\mathcal{I}' \le 0$) corresponds to the case where $x_{2,i}$ and $x_1$ are moving away (respectively $x_{2,i+1}$ and $x_1$).

Different cases can be mentioned:

- If $\mathcal{I}, \mathcal{I}' \le 0$, $x_{2,i}$ and $x_{2,i+1}$ are moving away from $x_1$, in consequence $f_2$ is moving away from $x_1$, no such singular contact can appear during the iteration.

- If $\mathcal{I} > 0$ and $\mathcal{I}' \le 0$: from $a = 0$, we deduce

$$ 
\begin{align}
& \Rightarrow (\overline{x_{2,i}x_1} + \overline{x_1x_{2,i+1}}) \wedge \overline{(v_{2,i+1} - v_1 - v_{2,i} + v_1)} = 0 \nonumber \\
& \Rightarrow (\overline{x_{2,i}x_1} + \overline{x_1x_{2,i+1}}) \wedge (\mathcal{I}'\overline{x_1x_{2,i+1}} - \mathcal{I}\overline{x_1x_{2,i}}) = 0 \nonumber \\
& \Rightarrow \mathcal{I}'(\overline{x_{2,i}x_1} \wedge \overline{x_1x_{2,i+1}}) - \mathcal{I}(\overline{x_1x_{2,i+1}} \wedge \overline{x_1x_{2,i}}) = 0 \tag{18} \\
& \Rightarrow \overline{x_1x_{2,i+1}} \wedge (-\mathcal{I}'\overline{x_2,i}x_1 - \mathcal{I}\overline{x_1x_{2,i}}) = 0 \nonumber \\
& \Rightarrow \overline{x_1x_{2,i+1}} \wedge \overline{x_1x_{2,i}}(\mathcal{I}' - \mathcal{I}) = 0. \nonumber
\end{align} $$
The colinearity condition between $\overline{x_1x_{2,i}}$ and $\overline{x_1x_{2,i+1}}$ implies also the colinearity with $\overline{x_1x_\alpha}$; and a colinearity between $(\overline{v_{2,i}} - \overline{v_1})$ and $(\overline{v_{2,i+1}} - \overline{v_1})$. Then by decomposition of (17), we have:

$$ t^*(\alpha) = \frac{\alpha ||\overline{x_1x_{2,i}}|| + (1-\alpha) ||\overline{x_1x_{2,i+1}}||}{\alpha ||\overline{(v_{2,i}-v_1)}|| + (1-\alpha) ||\overline{(v_{2,i+1}-v_1)}||} \left\langle \frac{\overline{x_1x_\alpha}}{||\overline{x_1x_\alpha}||}, \frac{\overline{(v_\alpha - v_1)}}{||\overline{(v_\alpha - v_1)}||} \right\rangle. \quad (19) $$

The shape of $\alpha \mapsto t^*(\alpha)$ is therefore as in Fig. 11.

<figure><img src="image_12.png" /><figcaption>Figure 11: Shape of $t^*(\alpha)$ if $\mathcal{I} > 0$ and $\mathcal{I}' \le 0$.</figcaption></figure>

As a consequence, $\min_{\alpha \in [0,1]}(t^*(\alpha) \ge 0) = t^*(\alpha = 1)$.

- Equivalently, if $\mathcal{I}' > 0$ and $\mathcal{I} \le 0$, $\min_{\alpha \in [0,1]}(t^*(\alpha) \ge 0) = t^*(\alpha = 0)$.

- Following the same strategy and decomposition, if $\mathcal{I}', \mathcal{I} > 0$, $\min_{\alpha \in [0,1]}(t^*(\alpha) \ge 0) = \min(\{t^*(\alpha = 0) \ge 0, t^*(\alpha = 1) \ge 0\})$.

- Moreover, if $\mathcal{I} = \mathcal{I}' \ne 0$, from (16), $c = b = 0$ implies

$$ \begin{aligned} \mathcal{I}(\overline{x_{2,i}x_1} + \overline{x_1x_{2,i+1}}) &= (\overline{v_1 - v_{2,i}}) + (\overline{v_{2,i+1} - v_1}) \\ \text{or } \mathcal{I}(\overline{x_{2,i}x_{2,i+1}}) &= (\overline{v_{2,i+1} - v_{2,i}}). \end{aligned} $$

A contact will appear between the three nodes $x_{2,i}$, $x_1$ and $x_{2,i+1}$ (i.e. the triangle shrinks to a point) at $t^*(\alpha) = \frac{\langle x_1(0) - x_{\alpha^*}(0), x_1(0) - x_{\alpha^*}(0) \rangle}{\mathcal{I}(x_1(0) - x_{\alpha^*}(0), x_1(0) - x_{\alpha^*}(0))} = \frac{1}{\mathcal{I}}$. But one must be recalled that contact is solved together with a balance law of the form (3): the CFL condition imposes a limited deformation of the cells during the iteration. Therefore $1/\mathcal{I} > \Delta t_{CFL}$ and no contact shall appears during the iteration (under this condition, we also respect $a = 0$ through the contact definition between $x_{2,i}$ and $x_{2,i+1}$).

The numerical strategy applied in §4 for the 2D case is described in Algorithm 1. This algorithm is designed for an arbitrary number of materials $M$. When intra-material contact (i.e., contact between parts of the same material) cannot occur,
**Algorithm 1** Two-dimensional contact detection algorithm.

for $\beta \in \{1, \dots, M\}$ (where $M$ is the number of materials that may be in contact)
do

$\partial\omega_{\bar{\beta}}$ is the master surface, $\partial\omega_{\beta}$ is the slave surface and $\bar{\beta} = \{1, \dots, M\}$

for $f_{\bar{\beta},i} \in \partial\omega_{\bar{\beta}}$ do

$f_{\bar{\beta},i}$ is defined by two points $x_{\bar{\beta},i}$ and $x_{\bar{\beta},i+1}$

for $x_{\beta,j} \in \partial\omega_{\beta}$ do

Solve $\Gamma_{i,j}(\alpha) = \alpha^2a + \alpha b + c = 0$, where $a, b, c$ are defined in (12)

if $a = b = c = 0$ then

Evaluate contact time between $x_{\beta,j}$ and points $x_{\bar{\beta},i}$ and $x_{\bar{\beta},i+1}$, denoted as
$t_c^{i,j}$

else

$$ \Delta = b^2 - 4ac $$

if $\Delta \ge 0$ then

Compute $\alpha^*$ such that $\Gamma_{i,j}(\alpha^*) = 0$

if $0 \le \alpha^* \le 1$ then

Contact will occur between $x_{\beta,j}$ and $x_{\alpha^*}$

Contact time $t_c^{i,j}$ is obtained from (14)

end if

end if

end if

end for

end for

$t^* = \min_{i,j}(t_c^{i,j}, t_c^{i,j} \ge 0)$ as the minimum of all possible contact times $t_c^{i,j}$

If $t^* \in [0, \Delta t]$, contact occurs within the time step

end for
$\bar{\beta} = \{1, \dots, M\}$ can be simplified to $\bar{\beta} = \{1, \dots, M\} \setminus \beta$. Thus, our formulation naturally accommodates intra-material contact. In the following discussion, we analyze Algorithm 1 for the case of two potentially impacting materials.

Algorithm 1 takes as input the surfaces of the objects of interest, where the master surface is described by its faces and the slave surface by its points. The algorithm iterates over all points and faces of these surface boundaries and performs the following steps:

- First, solves $\Gamma(\alpha) = 0$: if for a given point-face pair, a solution $\alpha^*$ exists in $[0, 1]$, contact may occur between them.

- If such an $\alpha^* \in [0, 1]$ exists, computes the contact time $t_c^{i,j}$ between point $x_{1,j} \in \partial\omega_1$ and face $f_{2,i} \in \partial\omega_2$.

- After testing all combinations, determines the contact time $t^*$ between $\omega_1$ and $\omega_2$ by finding the minimum of all computed $t_c^{i,j}$.

- Contact occurs within the time step if $t^* \in [0, \Delta t]$.

The complexity of this 2D generalization of the 1D procedure described in §3.1 is $O(N_1 \times N_2)$. For $M$ materials, this complexity can increase significantly. However, as in the 1D case, we can leverage a priori knowledge about which parts of surfaces $\partial\omega_1$ and $\partial\omega_2$ are more likely to impact (or conversely, unlikely to impact) to reduce complexity by appropriately defining the set $\bar{\beta}$. Additionally, since the computations of contact times $(t_c^{i,j})_{i \in \{1, \dots, N_1\}, j \in \{1, \dots, N_2\}}$ are independent, they can be efficiently parallelized, with only a simple min reduction needed for computing $t^*$.

### 3.3. The three-dimensional extension

The three-dimensional extension follows a similar approach: we consider two surfaces $\partial\omega_1(t)$ and $\partial\omega_2(t)$, each defined by either their nodes or their faces, such that $\partial\omega_\beta = \cup_{i=1}^{N_\beta} x_{\beta,i} = \cup_{i=1}^{N_\beta} f_{\beta,i}$ for $\beta \in \{1, 2\}$. Given a node $x_{1,j}(t) \in \partial\omega_1(t)$ and a face $f_{2,i'}(t) \in \partial\omega_2(t)$, we aim to detect potential contact between $x_{1,j}$ and $f_{2,i'}$. An example is shown in Figure 12.

<figure><img src="image_13.png" /><figcaption>Figure 12: Three dimensional problem: materials 1 and 2 are defined via their nodes $(x_{1,j})_{j \in \{1, \dots, N\}}$ for $\partial\omega_1$ and $(x_{2,j})_{j \in \{1, \dots, N'\}}$ for $\partial\omega_2$, $N, N' \in \mathbb{N}$. For a material $\beta \in \{1, 2\}$, the faces between four nodes $x_{\beta,i_1}, x_{\beta,i_2}, x_{\beta,i_3}, x_{\beta,i_4}$ is denoted by $f_{i'}^\beta$.</figcaption></figure>

Consider face $f_{2,i'}$ as a triangle (if $f_{2,i'}$ is a quadrilateral as shown in Figure 12, the procedure first decomposes it into triangles). The vertices of $f_{2,i'}$ are denoted by $x_{2,i_1}$,
$x_{2,i_2}$, and $x_{2,i_3}$. Each node in this system has its own velocity field: $v_{1,j}$, $v_{2,i_1}$, $v_{2,i_2}$, and $v_{2,i_3}$, respectively. This local configuration is illustrated in Figure 13.

<figure><img src="image_14.png" /><figcaption>Figure 13: The node $x_{1,j} \in \partial\omega_1$ and the surface $f_{2,i'} \in \partial\omega_2$.</figcaption></figure>

Following the approach from the previous section, we need to detect contact between $x_{1,j}$ and a point on face $f_{2,i'}$. To achieve this, we parameterize $f_{2,i'}$. For notational simplicity, in the following calculations, we denote $x_{1,j} = x$, $v_{1,j} = v$, and $f_{2,i'} = f$. The contact detection problem between node $x$ and surface $f$ becomes: find $x_\alpha = \sum_i \alpha_i x_i$, with $\alpha = (\alpha_1, \alpha_2, \alpha_3) \in [0, 1]^3$ and $\sum_i \alpha_i = 1$ such that

$$x_1(0) - x_\alpha(0) = \int_0^{t^*} [v_\alpha(s) - v_1(s)] \, ds.$$

Assuming constant velocities during the time iteration, as in the one- and two-dimensional cases, this equation becomes

$$x_1(0) - x_\alpha(0) = t^* [v_\alpha(0) - v_1(0)]. \quad (20)$$

This is equivalent to the collinearity condition:

$$\overrightarrow{x_1(0)x_\alpha(0)} \wedge (\overrightarrow{x_1(0)} + t^*v_1(0))(\overrightarrow{x_\alpha(0)} + t^*v_\alpha(0)) = 0.$$

The problem reduces to finding $\alpha$ such that $\Gamma(\alpha) = 0$, where $\Gamma \in \mathbb{R}^3 \times \mathbb{P}_2$ is defined analogously to (11) from the two-dimensional case by$^6$

$$\begin{aligned} \Gamma(\alpha) &= \frac{1}{t^*} \overrightarrow{x_1 x_\alpha} \wedge (\overrightarrow{x_1} + t^* v_1)(\overrightarrow{x_\alpha} + t^* v_\alpha), \\ &= \frac{1}{t^*} (\overrightarrow{x_1 x_\alpha} \wedge \overrightarrow{x_1 x_\alpha} + t^* \overrightarrow{x_1 x_\alpha} \wedge \overrightarrow{v_1 v_\alpha}) = \overrightarrow{x_1 x_\alpha} \wedge \overrightarrow{v_1 v_\alpha}, \end{aligned} \quad (21)$$

where $t^*$ is the contact time between $x_1$ and $x_\alpha$. Using the linear combinations $x_\alpha$ and $v_\alpha$ (since $v_\alpha = \sum_i \alpha_i v_{f_{2,i}}$) and noting that $\alpha_3 = 1 - \alpha_1 - \alpha_2$, we obtain

$$\begin{aligned} \Gamma(\alpha) &= \left( \sum_i \alpha_i \overrightarrow{x_1 x_{f_{2,i}}} \right) \wedge \left( \sum_i \alpha_i \overrightarrow{v_1 v_{f_{2,i}}} \right), \\ &= \alpha_1^2 (\overrightarrow{x_{f_{2,3}} x_{f_{2,i}}} \wedge \overrightarrow{v_{f_{2,3}} v_{f_{2,i}}}) + \alpha_2^2 (\overrightarrow{x_{f_{2,3}} x_{f_{2,2}}} \wedge \overrightarrow{v_{f_{2,3}} v_{f_{2,2}}}) \\ &\quad + \alpha_1 \alpha_2 (\overrightarrow{x_{f_{2,3}} x_{f_{2,1}}} \wedge \overrightarrow{v_{f_{2,3}} v_{f_{2,1}}}) + \overrightarrow{x_{f_{2,3}} x_{f_{2,2}}} \wedge \overrightarrow{v_{f_{2,3}} v_{f_{2,1}}}) \\ &\quad + \alpha_1 (\overrightarrow{x_1 x_{f_{2,3}}} \wedge \overrightarrow{v_{f_{2,3}} v_{f_{2,1}}}) + \overrightarrow{x_{f_{2,3}} x_{f_{2,1}}} \wedge \overrightarrow{v_1 v_{f_{2,3}}}) \\ &\quad + \alpha_2 (\overrightarrow{x_1 x_{f_{2,3}}} \wedge \overrightarrow{v_{f_{2,3}} v_{f_{2,2}}}) + \overrightarrow{x_{f_{2,3}} x_{f_{2,2}}} \wedge \overrightarrow{v_1 v_{f_{2,3}}}) \\ &\quad + (\overrightarrow{x_1 x_{f_{2,3}}} \wedge \overrightarrow{v_1 v_{f_{2,3}}}). \end{aligned} \quad (22)$$

<footnote>$^6$we omit the time dependence (0) for conciseness</footnote>
Assuming $\alpha^*$ exists and is unique such that $\Gamma(\alpha^*) = 0$, we can determine $x_{\alpha^*}$ and $t_c^*$. If $t_c^* < \Delta t$, contact will occur during the next time step.

Unlike the 2D case, $\alpha_1, \alpha_2 \mapsto \Gamma(\alpha_1, \alpha_2) = (\gamma_x(\alpha_1, \alpha_2), \gamma_y(\alpha_1, \alpha_2), \gamma_z(\alpha_1, \alpha_2))^T$ is a vector of bivariate quadratic functions. In 2D, while $\Gamma(\alpha)$ was also a vector, solving $\Gamma(\alpha) = 0$ reduced to solving $\gamma_z(\alpha) = 0$ (see (11)). Finding solutions to $\Gamma(\alpha) = 0$ in 3D is more complex. The procedure involves finding zeros of $\gamma_x(\alpha)$, $\gamma_y(\alpha)$, and $\gamma_z(\alpha)$, which correspond to the scalar product of $\Gamma(\alpha)$ projected onto the $\mathrm{Vect}(y, z)$, $\mathrm{Vect}(z, x)$, and $\mathrm{Vect}(x, y)$ planes, respectively. Contact may occur if and only if these functions share common zeros. In such cases, we compute the contact time and verify if it falls within the interval $[0, \Delta t]$.

In what follows, we detail this procedure step by step. Note that, unlike the 1D and 2D cases, we propose using an approximation of the contact time (and thus the contact position) in 3D. This choice is motivated by the non-negligible complexity of exactly solving the bivariate quadratic system $\Gamma(\alpha) = 0$. As discussed in Section 1 and illustrated in Figure 3, we only need first-order convergence with respect to $\Delta x$ for the contact time and position to ensure the expected convergence rate for the physical quantities of interest. Therefore, in the following section, we ensure first-order accuracy in the computation of $\alpha^*$ and $t^*$ to maintain the overall accuracy of the numerical method.

*Review of Bivariate Quadratic Functions.* A bivariate quadratic function takes the form

$$f(x, y) = Ax^2 + Bxy + Cy^2 + Dx + Ey + F, \forall (x, y) \in \mathbb{R}^2. \quad (23)$$

The solutions to $f(x, y) = 0$ define a conic section, which is the intersection of a surface with a plane. The nature of this intersection depends on the coefficients of $f$:

* if $B^2 - 4AC < 0$, the conic section is an ellipse (or a point in the degenerate case),

* if $B^2 - 4AC = 0$, the conic section is a parabola (or tangent lines in the degenerate case),

* if $B^2 - 4AC > 0$, the conic section is a hyperbola (or parallel lines in the degenerate case).

These cases are illustrated in Figure 14. For an ellipse, the extremum value (minimum or maximum, depending on the signs of $A$ and $C$) of the elliptic paraboloid occurs at

$$x_m = - \frac{2CD - EB}{4AC - B^2}, \quad (24)$$

$$y_m = - \frac{2AE - DB}{4AC - B^2}. \quad (25)$$

*Application to Contact Detection.* For our problem, we have

$$\gamma_p(\alpha_1, \alpha_2) = A\alpha_1^2 + B\alpha_1\alpha_2 + C\alpha_2^2 + D\alpha_1 + E\alpha_2 + F, \forall (\alpha_1, \alpha_2) \in \mathbb{R}^2, p \in \{x, y, z\}. \quad (26)$$

Note that $f_2 \in \partial\Omega$ represents a triangular closed and compact space. We distinguish two cases:
<figure><img src="image_15.png" /><figcaption>Figure 14: Types of non-degenerate conic sections: 1) circle, 2) ellipse, 3) parabola, 4) hyperbola.
(Source: Wikipedia)</figcaption></figure>

* When $B^2 - 4AC \ge 0$, the curve $\gamma_p(\alpha_1, \alpha_2) = 0$ represents a parabola, hyperbola, or lines. Consequently, if there exists $(\alpha_1, \alpha_2) \in f_2$ such that $\gamma_p(\alpha_1, \alpha_2) = 0$, then there must also exist $(\alpha_1, \alpha_2) \in \partial f_2$ satisfying $\gamma_p(\alpha_1, \alpha_2) = 0$. Therefore, we can evaluate zeros along the boundaries of $f_2$, effectively reducing this to a two-dimensional problem. In practice, we approximate the zeros of $\Gamma$ in $f_2$ by computing intersections of linear approximations of these conics on the surface $f_2$ (see the $\gamma_z = 0$ and $\gamma_x = 0$ lines in Figure 15). This ensures at least first-order accuracy in the approximation of the contact time and position.

* When $B^2 - 4AC < 0$, the curve $\gamma_p(\alpha_1, \alpha_2) = 0$ represents an ellipse or a point. In the case of an ellipse, every point $x$ on this ellipse satisfies the contact condition in the projected space—that is, every such point $x$ converges toward $x_1$ (the ellipse shrinks to a point, specifically $x_1$). We propose a first-order approximation of the contact time and position by approximating the zeros of $\Gamma$ in $f_2$ using the barycenter of the ellipse, computed from the extremum given by (25) (see the $\gamma_y = 0$ ellipse in Figure 15). While other approaches may exist, our numerical results in Section 4 demonstrate that this approximation is sufficient for consistency.
<figure><img src="image_16.png" /><figcaption>Figure 15: a) Distribution of exact solution sets of $\gamma_p(\alpha) = 0$ for $p \in \{x, y, z\}$; b) Distribution of solution sets obtained by the first-order approximation of $\gamma_p(\alpha) = 0$ for $p \in \{x, y, z\}$.</figcaption></figure>

Following the approach presented in Section 3.2 (14), the computation of $t^*$ is given by

$$t^* = t(\alpha^*) = \frac{\|\overrightarrow{x_1x_{\alpha^*}}\|}{\|(v_{\alpha^*} - v_1)\|} \left\langle \frac{\overrightarrow{x_1x_{\alpha^*}}}{\|(v_{\alpha^*} - v_1)\|}, \frac{\overrightarrow{(v_{\alpha^*} - v_1)}}{\|(v_{\alpha^*} - v_1)\|} \right\rangle, \quad (27)$$

$$= \frac{\langle x_1(0) - x_{\alpha^*}(0), x_1(0) - x_{\alpha^*}(0) \rangle}{\langle v_{\alpha^*}(0) - v_1(0), x_1(0) - x_{\alpha^*}(0) \rangle}. \quad (28)$$

**Note:** Equation (28) holds for any $N$-dimensional Euclidean space equipped with the $L_2$-norm (denoted by $\| \cdot \|$), where $N \in \mathbb{N}^*$.

The proposed 3D algorithm is detailed in Algorithm 2. As in the previous section on two-dimensional problems, Algorithm 2 is designed for an arbitrary number of materials $M$. The methodology in three-dimensional space remains consistent: Algorithm 2 takes as input the surfaces of the objects of interest. The master object ($\omega_{\bar{\beta}}$) is described by its faces, while the slave object ($\omega_{\beta}$) is described by its points. The algorithm iterates over all points and surfaces of the two volume boundaries and:

- First solves $\Gamma(\alpha) = 0$: this is performed with first-order approximation in three projected spaces: $\gamma_z(\alpha) = 0$ in space $(\vec{x}, \vec{y})$, $\gamma_x(\alpha) = 0$ in space $(\vec{y}, \vec{z})$, and $\gamma_y(\alpha) = 0$ in space $(\vec{z}, \vec{x})$. If solution sets exist in all three projected spaces and lie within the master surface, the contact point $x_c = x_1(t^*) = x_1(t(\alpha^*))$ is determined as the intersection of contact spaces in the three projected spaces.

- If $\alpha^* = (\alpha_1^*, \alpha_2^*, \alpha_3^*)$ exists within the surface, the contact time $t_c^{i,j}$ between point $x_{\beta,j} \in \partial\omega_{\beta}$ and surface $f_{\bar{\beta},i} \in \partial\omega_{\bar{\beta}}$ is computed.

- After testing all combinations, the contact time $t^*$ between $\omega_{\beta}$ and $\omega_{\bar{\beta}}$ is determined by finding the minimum of all computed $t_c^{i,j}$ values.

- Contact occurs within the time step if $t^* \in [0, \Delta t]$.
**Algorithm 2** Three-dimensional contact detection algorithm.

for $\beta \in \{1, \dots, M\}$ (where $M$ is the number of materials that may be in contact)
do

$\partial\omega_{\bar{\beta}}$ is the master surface, $\partial\omega_{\beta}$ is the slave surface and $\bar{\beta} = \{1, \dots, M\}$

for $f_{\bar{\beta},i} \in \partial\omega_{\bar{\beta}}$ do

$f_{\bar{\beta},i}$ is defined by three points $(x_{f_{\bar{\beta},1}}, x_{f_{\bar{\beta},2}}, x_{f_{\bar{\beta},3}})$

for $x_{\beta,j} \in \partial\omega_{\beta}$ do

for each component $\gamma_p$ of $\Gamma$, $p \in \{x, y, z\}$ do

Solve $\gamma_p(\alpha) = A\alpha_1^2 + B\alpha_1\alpha_2 + C\alpha_2^2 + D\alpha_1 + E\alpha_2 + F$ where $A, B, C, D, E, F$ are defined in (22)

$$ \Delta = B^2 - 4AC $$

if $\Delta < 0$ then

Compute $(x_m, y_m)$ from (25) (case of ellipse or point)

$$ \alpha_p = (\alpha_{p,1}, \alpha_{p,2}, \alpha_{p,3}) \leftarrow (x_m, y_m) $$

else if $\Delta \ge 0$ then

Evaluate zeros on $\partial f_{\bar{\beta},i}$ (case of parabola or hyperbola)

If zeros exist on $\partial f_{\bar{\beta},i}$: perform linear approximation of zeros on $f_{\bar{\beta},i}$

end if

end for

($\alpha_{c,1}, \alpha_{c,2}, \alpha_{c,3}$) $\leftarrow$ contact point, computed as the intersection of contact spaces in the three projected spaces

$t_c^{i,j} \leftarrow$ corresponding contact time from (28)

end for

$t^* = \min_{i,j}(t_c^{i,j}, t_c^{i,j} \le 0)$ as the minimum of all possible contact times $t_c^{i,j}$

If $t^* \in [0, \Delta t]$, contact occurs within the time step

end for

end for
The complexity of the 3D generalization of the 1D procedure described in Section 3.1 is $O(N_1 \times N_2)$, where $N_1$ and $N_2$ are the numbers of points on the surfaces of materials 1 and 2. While this complexity is the same as in 2D, the values of $N_1$ and $N_2$ can grow significantly faster in 3D as meshes are refined. Similarly, with $M$ potentially impacting materials, the computational cost can increase substantially. However, as in the 1D and 2D cases, the complexity can be reduced by leveraging a priori knowledge about which parts of surfaces $\partial\omega_1$ and $\partial\omega_2$ are more likely to come into contact. Furthermore, the computation of contact times $(t_c^{i,j})_{i\in\{1,...,N_1\},j\in\{1,...,N_2\}}$ remains independent, as in the 1D and 2D cases. This independence enables efficient parallelization in practice, with the computation of $t^*$ requiring only a simple min reduction operation.

# 4. Numerical tests

In this section, we evaluate the efficiency of our method through various multi-dimensional tests in $\mathbb{R}^n$ for $n \in \{1, 2, 3\}$. We begin in Section 4.1 by presenting the classical Node-to-Segment (NTS) detection algorithm from the literature. This a posteriori detection method has significant limitations that can severely impact convergence, as demonstrated in our benchmarks. Understanding these limitations provides context for comparing the performance of our proposed CETT approach. The NTS algorithm serves as our baseline for comparison. The implementation of the CETT approach used in these numerical tests is available in potatoes^7.

For our benchmarks, we focus on interactions between two materials, $\omega_1$ and $\omega_2$. While in some cases one material remains static while the other moves, we generally consider scenarios where both materials are in motion. For simplicity, we examine basic advection problems where, within a time step $[0, T]$, the velocity field of each node in both materials is known and constant. In Section 4.2, we demonstrate how the CETT approach inherently handles the pathological cases identified in Section 4.1. The subsequent sections provide quantitative analyses, including detailed convergence and performance studies. Section 4.3 examines two-dimensional benchmarks, while Section 4.4 focuses on three-dimensional cases. Additionally, we investigate the robustness of the CETT method across different mesh configurations in Section 4.4.2.

## 4.1. Review of the NTS Method for Two-Dimensional Problems [18]

In this section, we review the NTS method that serves as our baseline for comparison. We particularly focus on the limitations of this a posteriori detection strategy.

The NTS technique, introduced in Section 1, is an a posteriori approach that detects and corrects inter-penetrations that should not have physically occurred, often with limited information about the pre-contact dynamics. Consider a master surface $\partial\omega_2$ described by its faces and a slave surface $\partial\omega_1$ described by its nodes. The algorithm first identifies, for each slave node in $\partial\omega_1$, the closest node $x_{2,C}$ on the master surface $\partial\omega_2$. It is important to note that the NTS strategy pairs these closest nodes based solely on proximity, independent of the system dynamics. Once the closest node $x_{2,C}$ to $x_1 \in \partial\omega_1$ is found on $\partial\omega_2$, we consider it to lie on face $f_C^2$ defined by nodes $x_{2,L}$ and $x_{2,R}$.

<footnote n="7">⁷potatoes - Potatoid hydrOdynamics And conTAct deTechtiOn procedureES. CEA code developed by Sixtine Michel</footnote>
Specifically, $f_C^2$ is selected because it contains the closest point $x_{2,C} \in f_C^2 = [x_{2,L}, x_{2,R}]$. This configuration is illustrated in Figure 16. The NTS method then computes $g_R$ and $g_L$ as:

$$g_R = \langle \mathbf{g}, \overrightarrow{x_{2,C}x_{2,R}} \rangle, \qquad (29)$$

$$g_L = \langle \mathbf{g}, \overrightarrow{x_{2,C}x_{2,L}} \rangle, \qquad (30)$$

where $\mathbf{g}$ denotes the vector $\overrightarrow{x_{2,C}x_1}$.

<figure><img src="image_17.png" /><figcaption>Figure 16: Configuration of a slave node $x_1$ and the master surface $\partial\omega_2$.</figcaption></figure>

The quantities $g_L$ and $g_R$ represent the projections of $\mathbf{g}$ onto the face $f_C^2$ in two different orientations. Following the notation of Zavarise et al. [18], three cases can occur:

* $g_R \ge 0$ and $g_L < 0$: inside-right configuration

* $g_L \ge 0$ and $g_R < 0$: inside-left configuration

* $g_L \ge 0$ and $g_R \ge 0$: in-of-both configuration

Let $\mathbf{n}$ denote the outward normal of the relevant master surface (either $\overrightarrow{x_{2,C}x_{2,R}}$ or $\overrightarrow{x_{2,C}x_{2,L}}$). If $g_N = \langle \mathbf{g}, \mathbf{n} \rangle > 0$, we have an open-gap configuration with no penetration detected. If $g_N = \langle \mathbf{g}, \mathbf{n} \rangle \le 0$, we have a closed-gap configuration, indicating penetration or contact between the slave node and the master surface. We define $\mathbf{g}_N$ as the projection of $\mathbf{g}$ onto the outward normal $\mathbf{n}$: $\mathbf{g}_N = \langle \mathbf{g}, \mathbf{n} \rangle \mathbf{n}$.

<figure><img src="image_18.png" /><figcaption>Figure 17: Penetration detection between surfaces showing: a) inside-right configuration and b) in-of-both configuration. The slave node $x_{1,S}$ and master surface points $x_2$ are shown with their respective normal vectors.</figcaption></figure>

The NTS technique exhibits several pathological cases where its reliability becomes questionable. Since this method detects penetrations without considering the historical motion (i.e., the velocity fields during the time step), it cannot accurately determine either the exact impact point or the time of contact. To illustrate a practical and
potentially problematic case, consider the sharp-edged shape $\omega_1(t)$ shown on the left side of Figure 18. During the time step $[t, t + \Delta t]$, shape $\omega_1(t)$ moves from left to right and penetrates the static shape $\omega_2$. In this particular situation, the NTS method incorrectly identifies the entry point on the top surface because, at time $t + \Delta t$, node $x_1$ is closer to $x_{2,C}$ than to any point on the vertical surface of $\omega_2$, despite the actual penetration occurring through the vertical surface. While reducing the time step can

<figure><img src="image_19.png" /><figcaption>Figure 18: In this situation, the entering point is wrongly evaluated with NTS.</figcaption></figure>

mitigate the pathological situation shown in Figure 18, as node $x_1(t + \Delta t)$ will be positioned closer to the vertical face of $\omega_2$, allowing the NTS strategy to correctly identify the entry side. Such time step reduction naturally occurs when refining the mesh due to the CFL condition. However, if the refinement of $\omega_2$ proceeds faster than the decrease in $\Delta t$, the horizontal face will persistently be identified as the entry face, indicating that convergence to the correct solution cannot be guaranteed in certain scenarios.

A different problematic configuration is illustrated in Figure 19.

<figure><img src="image_20.png" /><figcaption>Figure 19: Case where the contact is not detected.</figcaption></figure>

In this case, penetration should be detected but is missed because material $\omega_1$ does not penetrate $\omega_2$ at either time $t$ or $t + \Delta t$, even though contact occurs at some intermediate time $t^* \in [t, t + \Delta t]$. While mesh refinement and the associated time step reduction (due to the CFL condition) will eventually lead to penetration detection at some critical time step, this requires unnecessarily intensive computational resources and may still not guarantee correct identification of the contact surface.

**Remark:** Although this paper focuses on contact detection, the subsequent correction phase merits attention. Currently, penetration correction is often performed without consideration of the system dynamics: once penetration is detected, materials are arbitrarily repositioned into contact along the normal between the penetrating node and the nearest surface. Our proposed CETT algorithm can complement existing NTS implementations by providing more accurate penetration correction. Since CETT precisely determines the contact location, it can compute the actual penetration direction (between the penetrating node and the contact point), enabling more physically accurate corrections than the conventional normal-based approach.

While these test cases may appear simplistic, they highlight important concerns in impact-contact calculations where such pathological situations may arise dynamically
during computation. The numerical tests presented in the following section demonstrate the robustness of the CETT technique through convergence analysis, particularly focusing on its effectiveness in handling the aforementioned pathological cases.

## 4.2. Application of the CETT Detection Strategy to the NTS Pathological Case

In this section, we examine how the CETT method handles the pathological case previously identified for NTS methods. We consider a simple test case involving two identical triangles $ω_1$ and $ω_2$, as illustrated in Figure 20. The initial conditions are defined as follows:

* $ω_1$: At $t = 0$, triangle $ω_1$ is defined by vertices $x_{1,1}(0) \cup x_{1,2}(0) \cup x_{1,3}(0) = (5.00, 0.00) \cup (2.00, 0.00) \cup (0.00, 0.00)$. Triangle $ω_1$ remains stationary throughout the computation, with zero velocity at all nodes. The material has constant density and no initial constraints.

* $ω_2$: At $t = 0$, triangle $ω_2$ is defined by vertices $x_{2,1}(0) \cup x_{2,2}(0) \cup x_{2,3}(0) = (4.90, 0.10) \cup (9.90, 0.10) \cup (9.90, 2.00)$. Triangle $ω_2$ moves with a uniform velocity $v_{2,i}(t_0) = (0, -1)$ for all nodes $i \in \{1, 2, 3\}$. The material has constant density and no initial constraints.

This benchmark allows for analytical computation of the contact time $t^* = 0.96$, based on the distance between $x_{2,2}$ and $ω_1$ given by $d = 0.1(1 - \frac{2}{5})$. The CETT algorithm accurately determines this contact time, regardless of mesh discretization. In contrast, the NTS strategy fails to detect contact unless an adequately small time step is used.

While the time step is typically governed by CFL conditions based on material sound velocities in $ω_1$ and $ω_2$, additional constraints may arise from contact detection requirements. For the NTS method to detect contact, it requires a time step smaller than $\Delta t_c \le t^*$, which implies a mesh size constraint of $\Delta x \le \|v_2\| \times \Delta t_c$. Although our example with $\Delta x = 1$ remains computationally feasible, scenarios involving higher velocities for $ω_2$ would impose prohibitively restrictive mesh size requirements.

<figure><img src="image_21.png" /><figcaption>Figure 20: Contact test case comparing NTS and CETT strategies. While NTS may miss the contact, CETT successfully detects it.</figcaption></figure>

It is important to note that even satisfying $\Delta t \le \Delta t_c$ may not guarantee proper contact detection. Consider a time step $\Delta t = \Delta t_c - \epsilon \le \Delta t_c$ where material $ω_2$ approaches but has not yet penetrated $ω_1$ at time $t = \Delta t_c$. If $\Delta t_c$ is not sufficiently small, $ω_2$ might pass through $ω_1$ without detection, resulting in the configuration shown in Figure 20 at time $t = 2\Delta t_c$.
<figure><img src="image_22.png" /><figcaption>Figure 21: Unstructured quadrilateral mesh of a circle with an average cell characteristic length of $\Delta x_{\text{mean}} = 0.16$.</figcaption></figure>

For this particular test case, CETT achieves exact contact time and location determination because all surfaces of $\omega_1$ and $\omega_2$ are planar. In subsequent sections, we examine spherical geometries where mesh refinement affects the accuracy of both methods.

### 4.3. Two-Dimensional Numerical Tests

This section presents our two-dimensional numerical tests. In Section 4.3.1, we examine a test case involving the impact between two circles. With circular geometries, the CETT strategy no longer provides exact solutions (as increasingly fine meshes are required to accurately approximate the circles), and mesh resolution becomes significant. This allows us to evaluate and compare the convergence properties and computational complexity of both CETT and NTS methods. In Section 4.3.2, we analyze a 2D configuration that reduces to a 1D case, verifying that CETT behaves as designed by naturally degenerating to the treatment described in Sections 1 and 3.1.

#### 4.3.1. Convergence and Complexity Studies: Circle-to-Circle Test Case

In this section, we analyze a test case involving contact between two circles, $\omega_1$ and $\omega_2$. This configuration is chosen because the tessellation of circles introduces numerical approximations (unlike the exact representation possible with polygons). We discretize the boundaries $\partial\omega_1$ and $\partial\omega_2$ using unstructured quadrilateral meshes, as shown in Figure 21.

This problem allows us to quantify the mesh's effect on the contact time $t^*$ and evaluate the convergence rates of both CETT and NTS methods with respect to this observable. We define $\omega_1(0)$ and $\omega_2(0)$ as circles of radius $r_1 = r_2 = 1$, with initial centers at $x_{r1} = (0,0)^T$ and $x_{r2} = (27,9)^T$ respectively. The velocity of each node within each material remains constant before contact, with $v_1(x,t) = v_1 = (3,1)^T$ for all $(x,t) \in \omega_1 \times [0,t^*]$ and $v_2(x,t) = v_2 = (-3,-1)^T$ for all $(x,t) \in \omega_2 \times [0,t^*]$. Figure 22 illustrates this test case configuration.
<figure><img src="image_23.png" /><figcaption>Figure 22: Configuration of the convergence test case described in Section 4.3.1.</figcaption></figure>

The initial positions and velocities were chosen to facilitate analytical determination of the contact time and position. The initial distance between the two contact points is $d_0 = \sqrt{27^2 + 9^2} - (r_1 + r_2)$. Given their relative velocity $v_1 - v_2$, the contact time is $t^* = \frac{d_0}{\|v_1 - v_2\|} \approx 4.183772234$. The contact position is then $x_1(t^*) = x_2(t^*) = (13.5, 4.5)^T$.

To evaluate the computational cost of both algorithms, we perform simulations using different mesh sizes with average characteristic lengths $\Delta x_{\text{mean}} \in \{0.32, 0.16, 0.08, 0.04, 0.02, 0.01\}$. We now compare the performance of the CETT technique with the NTS algorithm across these mesh resolutions.

<table><caption>Table 1: Comparison of CETT and NTS methods: convergence analysis and computational performance.</caption><thead><tr><th rowspan="2">Grid size<br>Δx</th><th colspan="3">CETT</th><th colspan="3">NTS</th></tr><tr><th>L<sub>1</sub> Error<br>on t*</th><th>Detection</th><th>CPU time</th><th>L<sub>1</sub> Error<br>on t*</th><th>Detection</th><th>CPU time</th></tr></thead><tbody><tr><td>0.16</td><td>0.000179714</td><td>Yes</td><td>0.09</td><td>N/A</td><td>Yes</td><td>0.08</td></tr><tr><td>0.08</td><td>8.51833e-05</td><td>Yes</td><td>0.35</td><td>N/A</td><td>Yes</td><td>0.24</td></tr><tr><td>0.04</td><td>3.80245e-05</td><td>Yes</td><td>0.96</td><td>N/A</td><td>Yes</td><td>1.23</td></tr><tr><td>0.02</td><td>1.44558e-05</td><td>Yes</td><td>6.33</td><td>N/A</td><td>Yes</td><td>4.99</td></tr><tr><td>0.01</td><td>3.61606e-06</td><td>Yes</td><td>42.09</td><td>N/A</td><td>Yes</td><td>22.8</td></tr><tr><td>Order</td><td>1.38</td><td>—</td><td>2.1</td><td>N/A</td><td>—</td><td>2.2</td></tr></tbody></table>

Table 1 compares the performance of both detection strategies in terms of accuracy and computational efficiency. While CETT provides quantifiable error measurements for contact time $t^*$, NTS only indicates whether contact is detected. Consequently, the convergence rate for the NTS method cannot be computed a posteriori.

The CETT method exhibits a convergence order of approximately 1.38, which lies between first and second order. This warrants further investigation: in our convergence study, we specified an average mesh size $\Delta x$ as input to the gmsh mesh generator [46]. This approach ensures relatively uniform meshes with comparable cell areas (see Figure 21), which is important for numerical simulations. However, it does not explicitly control the discretization of boundaries, which is crucial for accurate contact detection.

To illustrate this point, we conducted additional tests using an alternative meshing strategy: instead of using the average mesh size criterion, we divided the circle perimeters into 16, 32, 64, 80, 100, 120, 128, 140, and 256 segments. While the resulting
meshes are less optimal for numerical simulations due to the presence of high-aspect-ratio cells (see Figure 23, left), this boundary-focused discretization yields convergence rates much closer to second order. Table 2 presents these results, showing errors in contact time $t^*$ for different perimeter subdivisions. Figure 23 (right) demonstrates the observed second-order convergence rate for this alternative approach.

This second-order convergence is particularly significant because standard numerical methods in hydrodynamics typically achieve only first-order accuracy when shocks are present, indicating that CETT maintains the scheme's numerical accuracy.

<table><caption>Table 2: Convergence analysis of the CETT method using circle perimeter subdivision.</caption><thead><tr><th rowspan="2">Perimeter Division</th><th colspan="3">CETT</th></tr><tr><th><i>h</i></th><th><i>L</i><sub>1</sub> Error on <i>t</i><sup>*</sup></th><th>Detection</th></tr></thead><tbody><tr><td>16</td><td>0.3902</td><td>0.003621515821</td><td>Yes</td></tr><tr><td>32</td><td>0.1960</td><td>0.001406046652</td><td>Yes</td></tr><tr><td>64</td><td>0.0981</td><td>0.0003054220289</td><td>Yes</td></tr><tr><td>80</td><td>0.0713</td><td>0.0002014600199</td><td>Yes</td></tr><tr><td>100</td><td>0.0604</td><td>0.0001267324388</td><td>Yes</td></tr><tr><td>120</td><td>0.0491</td><td>9.410441495e-05</td><td>Yes</td></tr><tr><td>128</td><td>0.0462</td><td>1.161533753e-05</td><td>Yes</td></tr><tr><td>140</td><td>0.0436</td><td>7.04740544e-05</td><td>Yes</td></tr><tr><td>258</td><td>0.0238</td><td>2.235791962e-05</td><td>Yes</td></tr><tr><td>Order</td><td></td><td>2.01</td><td>—</td></tr></tbody></table>

The notable feature in the convergence curve (Figure 23, right) warrants explanation: the observed accuracy peak corresponds to a specific mesh configuration illustrated in Figure 24. This occurs when the mesh exhibits a particular symmetry where at time $t^*$, a slave node $x_{1,j} \in \partial\omega_1$ contacts a master surface $f_i^2$ of $\omega_2$, while simultaneously, a slave node $x_{2,i} \in \partial\omega_2$ contacts a master surface $f_j^1$ of $\omega_1$. In this symmetric configuration, the exact contact position $x_c = (13.5, 4.5)$ coincides (to machine precision) with the barycenter of these two slave nodes, leading to exceptionally accurate detection. This peak in accuracy occurs when the boundary nodes of the materials align closely with the true contact points on the spherical surfaces.

Regarding computational efficiency, CETT exhibits slightly higher CPU times than NTS. However, both methods share the same complexity order of approximately 2, indicating similar computational growth rates with mesh refinement. Thus, CETT achieves higher accuracy and better convergence properties while maintaining the same algorithmic complexity as the standard NTS method. Moreover, the additional computational cost of CETT remains modest and practically manageable.

### 4.3.2. Remark about 2D test-case degenerating toward a 1D test-case

In the previous section, we observed that CETT can exhibit distinct peaks in accuracy under certain conditions. Here, we further explore these accuracy characteristics by demonstrating that CETT, by design, achieves exact solutions for one-dimensional test cases. This property suggests that our algorithm successfully generalizes the one-dimensional approach described in the introduction.
<figure><img src="image_24.png" /><figcaption>Figure 23: Left: Quadrilateral mesh of a circle with boundary discretization of 64 segments. Right: Contact time error $t_c$ for varying boundary discretizations.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>Figure 24: Pattern of the contact at $t^*$ for the mesh generating the "dent"/"peak of accuracy" in the convergence study of Fig. 23-right. The (slave) node $x_{1,j} \in \partial\omega_1$ is detected in contact with the (master) surface $f_i^2$ composing $\omega_2$. Symetrically the (slave) node $x_{2,i} \in \partial\omega_2$ is detected in contact with the (master) surface $f_j^1$ composing $\omega_1$. The exact contact position $x_c = (13.5, 4.5)$ is the barycenter (machine precision) of these two nodes $x_{1,i}$ and $x_{2,j}$.</figcaption></figure>

To illustrate this point, Figure 25 presents a one-dimensional test configuration where $\omega_1 = [-1, 1]$ moves with a constant velocity $v_1 = -0.5$ along the $\vec{x}$ direction, while $\omega_2 = [5, 7]$ moves with a constant velocity $v_2 = -1$ in the same direction. The contact occurs at time $t_c = 8$ at position $x_c = (-3, 0)^T$. As expected, the *potatoes* code computes these values exactly.

<figure><img src="image_26.png" /><figcaption>Figure 25: One-dimensional problem.</figcaption></figure>

This exactness property extends to two-dimensional cases when dealing with planar boundaries, specifically when the discretized mesh $\omega_{1,h}$ perfectly matches the physical body $\omega_1$ (i.e., when $\omega_1 \setminus \omega_{1,h} = \emptyset$). In such cases, CETT provides exact contact detection. Figure 26 demonstrates this property through a two-dimensional extension of the previously discussed one-dimensional test case. Here, the contact occurs along the segment $[(-3, 0), (-3, 1)]$, representing the interaction between the right face of $\omega_1$ and the left face of $\omega_2$.
<figure><img src="image_27.png" /><figcaption>Figure 26: Two-dimensional contact test case with exact mesh-body correspondence.</figcaption></figure>

In this two-dimensional configuration, both the contact time and position are determined with perfect accuracy.

### 4.4. Three-Dimensional Numerical Tests

This section extends our analysis to three-dimensional problems. Following the structure of the previous section, we first evaluate CETT's convergence and computational complexity using a test case with analytically computable contact time and position (Section 4.4.1). Then, in Section 4.4.2, we assess the method's robustness across different quadrilateral boundary discretizations.

#### 4.4.1. Sphere-to-Sphere Contact: Convergence and Complexity Analysis

This test case represents a three-dimensional extension of the two-dimensional problem presented in Section 4.3.1, allowing us to evaluate CETT's convergence properties.

The system consists of two spheres, $ω_1$ and $ω_2$, each with radius $r_1 = r_2 = 1$ and initial centers at $x_{r1} = (0, 0, 0)^T$ and $x_{r2} = (27, 9, 9)^T$ respectively. Each point within $ω_1$ and $ω_2$ moves with constant velocities $v_1(x, t) = v_1 = (3, 1, 1)^T$ for all $(x, t) \in ω_1 \times \mathbb{R}_+$ and $v_2(x, t) = v_2 = (-3, -1, -1)^T$ for all $(x, t) \in ω_2 \times \mathbb{R}_+$. Figure 27 illustrates this configuration.

<figure><img src="image_28.png" /><figcaption>Figure 27: Three-dimensional contact test case configuration described in Section 4.4.1.</figcaption></figure>

Following calculations similar to those in Section 4.3.1, we determine the contact position as $x_c = (13.5, 4.5, 4.5)^T$ and the contact time as $t^* \approx 4.198488655$. The distance between either $∂ω_1$ or $∂ω_2$ and $x_c$ is $d_c = \sqrt{13.5^2 + 4.5^2 + 4.5^2 - 1}$, yielding $t^* = \frac{\sqrt{13.5^2+4.5^2+4.5^2-1}}{\sqrt{3^2+1^2+1^2}}$.

For mesh generation, we employ a hexahedral subdivision of the domains, resulting in quadrilateral boundary elements, as shown in Figure 28.
<figure><img src="image_29.png" /><figcaption>Figure 28: Unstructured hexahedral mesh of a sphere with average characteristic length $\Delta x_{\text{mean}} = 0.16$.</figcaption></figure>

Algorithm 2 in Section 3.3 describes the master surface using triangular faces. To accommodate our quadrilateral mesh elements, we must subdivide them into triangles. This subdivision is not unique (see Figure 29), and we will evaluate CETT's robustness with respect to different subdivision strategies.

<figure><img src="image_30.png" /><figcaption>Figure 29: Alternative triangular subdivisions of a quadrilateral element.</figcaption></figure>

For our initial convergence analysis, we employ subdivision scheme *a*) from Figure 29. We evaluate CETT's computational performance across different mesh resolutions, using average characteristic lengths $\Delta x_{\text{mean}} \in \{0.32, 0.16, 0.08, 0.04, 0.02, 0.01\}$.

<table><caption>Table 3: Convergence analysis of CETT using subdivision scheme *a*).</caption><thead><tr><th rowspan="2">Grid size<br>Δx</th><th colspan="2">CETT</th></tr><tr><th>L<sub>1</sub> Error</th><th>CPU time</th></tr></thead><tbody><tr><td>0.64</td><td>0.0016152</td><td>21.58</td></tr><tr><td>0.32</td><td>0.0006299210</td><td>138.28</td></tr><tr><td>0.16</td><td>0.00059724</td><td>1838.57</td></tr><tr><td>0.08</td><td>0.00014107</td><td>91301.25</td></tr><tr><td>Order</td><td>1.06</td><td>3.98</td></tr></tbody></table>

Table 3 presents the $L_1$ error in contact time $t^*$ and associated computational costs using subdivision scheme *a*). The method achieves first-order convergence, which aligns with the accuracy of the PDE numerical schemes discussed in Section 1. The computational complexity grows as $O(n^4)$, suggesting that CETT's complexity scales as $O(n^{d-1})$ in $d$ dimensions.

### 4.4.2. Sensitivity Analysis: Impact of Quadrilateral Subdivision

As discussed earlier, implementing Algorithm 2 with hexahedral meshes requires triangulating the boundary elements. Table 4 compares CETT's accuracy using the two subdivision schemes shown in Figure 29.
<table><caption>Table 4: Accuracy comparison between subdivision schemes a) and b) (Figure 29).</caption><thead><tr><th>Grid size<br>Δx</th><th>CETT - subdivision a)<br>l<sub>1</sub> Error</th><th>CETT - subdivision b)<br>l<sub>1</sub> Error</th></tr></thead><tbody><tr><td>0.64</td><td>0.0016152</td><td>0.00129991</td></tr><tr><td>0.32</td><td>0.0006299210</td><td>0.000540511</td></tr><tr><td>0.16</td><td>0.00059724</td><td>7.80684e-05</td></tr><tr><td>0.08</td><td>0.00014107</td><td>0.000164403</td></tr><tr><td>Order</td><td>1.06</td><td>1.17</td></tr></tbody></table>

Both subdivision schemes achieve similar convergence rates: 1.06 for scheme a) and 1.17 for scheme b). This three-dimensional convergence rate differs from the two-dimensional case (see Table 1), possibly due to the approximations introduced in representing two-dimensional cell surfaces. As in the two-dimensional case, the results in Table 4 show errors and convergence rates based on an average mesh size Δx specified to the mesh generator. While this approach produces regular meshes with uniform cell volumes suitable for numerical simulations, it does not guarantee optimal boundary discretization for contact detection. Consequently, the convergence rates in Table 4 may be conservative, as demonstrated in Section 4.3.1. This highlights the importance of mesh generation strategy: for contact-impact problems, we must balance boundary discretization quality with volumetric mesh properties.

Nevertheless, even with first-order convergence, the three-dimensional CETT extension maintains the numerical accuracy of the scheme, which is promising for contact/impact simulations. We note that subdivision scheme b) incurs higher computational costs, requiring approximately twice the number of operations compared to scheme a).

## 5. Conclusion

In this paper, we have presented a new contact detection algorithm called CETT (Continuous Extension of Thalès' Theorem). It is a 2D/3D extension of an efficient 1D contact detection algorithm described in §1. CETT provides a converging first/second order approximation of the contact time and location in 1D/2D/3D, with the order depending on the quality of the boundary meshes. In the numerical examples presented in this paper, the expected convergence rates were experimentally verified using original analytical benchmarks. The obtained convergence rates are compatible with those of classical numerical schemes used to simulate elasto-plastic materials in Lagrangian coordinates (which are generally of first order).

The new CETT method has been compared to the baseline approach called NTS (Node to Segment). Several limitations of the NTS method have been highlighted through simple pathological analytical cases where CETT performed successfully. Comparisons of the two methods for two-dimensional test cases show that the NTS technique and CETT are comparable in terms of CPU-time and complexity. While NTS remains slightly less computationally intensive, it provides significantly less information. The additional cost of CETT is justified by its benefits: it provides accurate information based on the dynamics and inherently avoids complex pathological cases that may arise during simulation.
An important feature of this new contact detection procedure is that it can be used
to enhance contact detection-correction methods implemented in legacy codes. For
instance, in a code where the NTS approach is implemented, CETT can be integrated
between the penetration detection and correction steps to evaluate the contact position
between the master surface and the slave node. Subsequently, the *non-penetration* cor-
rection can be applied in this contact direction instead of the (arbitrary) face outward
normal direction as traditionally done in the literature. We believe that CETT has the
potential to both improve existing algorithms and serve as a foundation for new ones.

Finally, the CETT algorithms in 1D, 2D, and 3D are all derived from the same
continuous procedure. These algorithms are inherently parallelizable and adaptable
for Lagrangian staggered and cell-centered schemes, while remaining largely indepen-
dent of the specific elasto-plastic model being used.

Acknowledgment

This work was performed within the Post doctoral project of Sixtine Michel at the
CEA CESTA.

Declarations

**Funding** Sixtine Michel, Gael Poette and Nicolas Therme were funded by the CEA CESTA.

**Conflicts of interest/Competing interests** The authors certify that there is no
actual or potential conflict of interest in relation to this article.
## References

[1] I. Peshkov, W. Boscheri, R. Loubère, E. Romenski, M. Dumbser, Theoretical and numerical comparison of hyperelastic and hypoelastic formulations for eulerian non-linear elastoplasticity, Journal of Computational Physics 387 (2019) 481–521.
doi:https://doi.org/10.1016/j.jcp.2019.02.039.
URL https://www.sciencedirect.com/science/article/pii/S0021999119301561

[2] G. Kluth, B. Després, Discretization of hyperelasticity on unstructured mesh with a cell-centered lagrangian scheme, Journal of Computational Physics 229 (24) (2010) 9092–9118.
doi:https://doi.org/10.1016/j.jcp.2010.08.024.
URL https://www.sciencedirect.com/science/article/pii/S0021999110004717

[3] M. L. Wilkins, Calculation of elastic-plastic flow. (1 1969).
URL https://www.osti.gov/biblio/4775245

[4] N. Favrie, S. Gavrilyuk, A Well-posed Hypoelastic Model Derived From a Hyperelastic One, John Wiley & Sons, Ltd, 2018, Ch. 14, pp. 417–427.
arXiv: https://onlinelibrary.wiley.com/doi/pdf/10.1002/9781119579311.ch14,
doi:https://doi.org/10.1002/9781119579311.ch14.
URL https://onlinelibrary.wiley.com/doi/abs/10.1002/9781119579311.ch14

[5] Dassault Systèmes, Abaqus 2022 Documentation, version 2022 (2022).
URL https://www.3ds.com/fr/produits-et-services/simulia/produits/abaqus/derniere-version/

[6] Commissariat à l'Énergie Atomique et aux Énergies Alternatives (CEA) et Centre Commun de Recherche de la Commission Européenne (JRC), Europlexus User's Manual, version 2024.0 (2024).
URL http://www-epx.cea.fr/index.php/epx-docs/25-epx-user-manual

[7] Ansys, Inc., LS-DYNA Keyword User's Manual, version R15.0 (2024).
URL https://www.dynasupport.com/manuals/ls-dyna-manuals

[8] Altair Engineering, Altair Radioss 2024 Documentation, version 2024 (2024).
URL https://2024.help.altair.com/2024/hwsolvers/rad/index.html

[9] R. D. Richtmyer, K. W. Morton, Difference methods for initial-value problems / Robert D. Richmyer, K.W. Morton., 2nd Edition, Interscience tracts in pure and applied mathematics ; 4, Interscience Publishers, New York, 1967.

[10] A. B. F. S. R. M., Methods in Computational Physics : Advances in Research and Applications, vol 3 Edition, 1964.

[11] D. Burton, Exact Conservation of Energy and Momentum in Staggered-grid Hydrodynamics With Arbitrary Connectivity, Vol. 395, 1990, pp. 7–19.
doi: 10.1007/3-540-54960-9_34.
[12] E. Caramana, D. Burton, C. Rousculp, M. Shashkov, P. Whalen, The construction of compatible hydrodynamics algorithms in multiple dimensions utilizing conservation of total energy, la-ur-99-2930, 1999. doi:10.2514/6.1999-3818.

[13] A. Bauer, D. Burton, E. Caramana, R. Loubère, M. Shashkov, P. Whalen, The internal consistency, stability, and accuracy of the discrete, compatible formulation of lagrangian hydrodynamics, Journal of Computational Physics 218 (2006) 572–593. doi:10.1016/j.jcp.2006.02.024.

[14] C. Mazeran, Sur la structure mathématique et l'approximation numérique de l'hydrodynamique lagrangienne bidimensionnelle, Ph.D. thesis, thèse de doctorat dirigée par Bruneau, Charles-Henri et Després, Bruno Informatique et mathématiques appliquées Bordeaux 1 2007 (2007). URL http://www.theses.fr/2007BOR13470

[15] R. Loubère, P.-H. Maire, B. Rebourcet, Chapter 13 - staggered and collocated finite volume schemes for lagrangian hydrodynamics, in: R. Abgrall, C.-W. Shu (Eds.), Handbook of Numerical Methods for Hyperbolic Problems, Vol. 17 of Handbook of Numerical Analysis, Elsevier, 2016, pp. 319–352. doi:https://doi.org/10.1016/bs.hna.2016.07.003.
URL https://www.sciencedirect.com/science/article/pii/S1570865916300059

[16] M. Gurtin, E. Fried, L. Anand, The Mechanics and Thermodynamics of Continua, 2010. doi:10.1017/CBO9780511762956.

[17] N. Burago, V. Kukudzhanov, A review of contact algorithms, Mechanics of Solids (2005) 44–85.

[18] G. Zavarise, L. De Lorenzis, The node-to-segment algorithm for 2d frictionless contact: Classical formulation and special cases, Computer Methods in Applied Mechanics and Engineering - COMPUT METHOD APPL MECH ENG 198 (2009) 3428–3451. doi:10.1016/j.cma.2009.06.022.

[19] G. Zavarise, L. De Lorenzis, The contact patch test for linear contact pressure distributions in 2d frictionless contact, Lecture Notes in Applied and Computational Mechanics 58 (01 2011). doi:10.1007/978-3-642-22167-5_5.

[20] G. Clair, B. Despres, E. Labourasse, A new method to introduce constraints in cell-centered lagrangian schemes, Computer Methods in Applied Mechanics and Engineering 261-262 (2013) 56–65. doi:https://doi.org/10.1016/j.cma.2013.03.013.
URL https://www.sciencedirect.com/science/article/pii/S0045782513000649

[21] J. A. González, K. Park, C. A. Felippa, R. Abascal, A formulation based on localized lagrange multipliers for bem–fem coupling in contact problems, Computer Methods in Applied Mechanics and Engineering 197 (6) (2008) 623–640. doi:https://doi.org/10.1016/j.cma.2007.08.020.
URL https://www.sciencedirect.com/science/article/pii/S0045782507003623

[22] J. Gardelle, B. Cassany, N. Szalek, A. Galtié, B. Bicrel, D. Hébert, I. Bertron, K. Pepitone, High power electron beam interaction with an aluminum target: Measurements and simulations, Journal of Applied Physics 126 (15) (2019) 153301. arXiv:https://pubs.aip.org/aip/jap/article-pdf/doi/10.1063/1.5123031/13019723/153301\_1\_online.pdf, doi:10.1063/1.5123031. URL https://doi.org/10.1063/1.5123031

[23] D. Hébert, L. Hallo, L. Voisin, T. Desanlis, A. Galtié, B. Bicrel, C. Mau-nier, P. Mercier, G. Duchateau, A KDP equation of state for laser-induced damage applications, Journal of Applied Physics 109 (12) (2011) 123527. arXiv:https://pubs.aip.org/aip/jap/article-pdf/doi/10.1063/1.3600210/13285144/123527\_1\_online.pdf, doi:10.1063/1.3600210. URL https://doi.org/10.1063/1.3600210

[24] D. Hébert, G. Seisson, J.-L. Rullier, I. Bertron, L. Hallo, J.-M. Chevalier, C. Thessieux, F. Guillet, M. Boustié, L. Berthe, Hypervelocity impacts into porous graphite: experiments and simulations, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 375 (2085) (2017) 20160171.

[25] B. Jodar, D. Hébert, B. Aubert, E. Andò, J.-L. Rullier, I. Bertron, Impacts into porous graphite: An investigation on crater formation and ejecta dis-tribution, International Journal of Impact Engineering 152 (2021) 103842. doi:https://doi.org/10.1016/j.ijimpeng.2021.103842. URL https://www.sciencedirect.com/science/article/pii/S0734743X21000294

[26] G. Seisson, D. Hébert, I. Bertron, J.-M. Chevalier, L. Hallo, E. Le-scoute, L. Videau, P. Combis, F. Guillet, M. Boustie, L. Berthe, Dy-namic cratering of graphite: Experimental results and simulations, In-ternational Journal of Impact Engineering 63 (2014) 18-28. doi:https://doi.org/10.1016/j.ijimpeng.2013.08.001. URL https://www.sciencedirect.com/science/article/pii/S0734743X13001565

[27] P. Pradel, F. Malaise, T. de Resseguier, C. Delhomme, G. Le Blanc, J. Ques-sada, Dynamic compaction of polyurethane foam: experiments and modelling, The European Physical Journal. Special Topics 227 (1-2) (2018) 3-16. doi:10.1140/epjst/e2018-00116-7. URL https://hal.science/hal-02301795

[28] R. J. LeVeque, Finite Volume Methods for Hyperbolic Problems, Cambridge Uni-versity Press, 2002. URL http://www.clawpack.org/fvmhp_materials/
[29] J. VonNeumann, R. D. Richtmyer, A Method for the Numerical Calculation of Hydrodynamic Shocks, Journal of Applied Physics 21 (3) (1950) 232-237. arXiv:https://pubs.aip.org/aip/jap/article-pdf/21/3/232/18309197/232_1_online.pdf, doi:10.1063/1.1699639. URL https://doi.org/10.1063/1.1699639

[30] N. Kikuchi, J. T. Oden, Contact Problems in Elasticity: A Study of Variational Inequalities and Finite Element Methods, SIAM Studies in Applied Mathematics, Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 1988. doi:10.1137/1.9781611970845.

[31] T. J. R. Hughes, R. L. Taylor, Unconditionally stable algorithms for structural dynamics analysis employing a mixed variational principle, Computer Methods in Applied Mechanics and Engineering 10 (2) (1978) 201-220.

[32] P. Alart, A. Curnier, A mixed formulation for frictional contact problems prone to newton like solution methods, Computer Methods in Applied Mechanics and Engineering 92 (3) (1991) 353-375.

[33] J. C. Simo, T. A. Laursen, An augmented lagrangian treatment of contact problems involving friction, Computer & Structures 42 (1) (1992) 97-116.

[34] F. D. Fischer, R. Rentzsch, On contact problems with rigid obstacles by penalty and lagrange-multiplier methods, International Journal for Numerical Methods in Engineering 40 (4) (1997) 667-681.

[35] M. Tur, F. Fuenmayor, P. Wriggers, A mortar-based frictional contact formulation for large deformations using lagrange multipliers, Computer Methods in Applied Mechanics and Engineering 198 (37) (2009) 2860-2873. doi:https://doi.org/10.1016/j.cma.2009.04.007. URL https://www.sciencedirect.com/science/article/pii/S0045782509001595

[36] R. Weyler, J. Oliver, T. Sain, J. Cante, On the contact domain method: A comparison of penalty and lagrange multiplier implementations, Computer Methods in Applied Mechanics and Engineering 205-208 (2012) 68-82, special Issue on Advances in Computational Methods in Contact Mechanics. doi:https://doi.org/10.1016/j.cma.2011.01.011. URL https://www.sciencedirect.com/science/article/pii/S0045782511000120

[37] G. Clair, B. Despres, E. Labourasse, A one-mesh method for the cell-centered discretization of sliding, Computer Methods in Applied Mechanics and Engineering 269 (2014) 315-333. doi:https://doi.org/10.1016/j.cma.2013.10.020. URL https://www.sciencedirect.com/science/article/pii/S0045782513002727

[38] D. Benson, J. Hallquist, A single surface contact algorithm for the post-buckling analysis of shell structures, Computer Methods in Applied Mechanics and Engineering 78 (2) (1990) 141-163. doi:10.1016/0045-7825(90)90072-8.
[39] L. Taylor, D. Flanagan, PRONTO 2D: A Two-Dimensional Transient Solid Dynamics Program, Tech. Rep. SAND86-0594, Sandia National Laboratories (1987).
doi:10.2172/6671798.
URL https://www.osti.gov/biblio/6671798

[40] Z.-H. Zhong, Finite Element Procedures for Contact-Impact Problems, Oxford Science Publications, Oxford University Press, Oxford, UK, 1993.

[41] A. B. Chaudhary, K.-J. Bathe, A solution method for static and dynamic analysis of three-dimensional contact problems with friction, Computers & Structures 24 (6) (1986) 855-873. doi:10.1016/0045-7949(86)90007-4.

[42] D. J. Benson, Computational methods in lagrangian and eulerian hydrocodes, Computer Methods in Applied Mechanics and Engineering 99 (2) (1992) 235-394.
doi:https://doi.org/10.1016/0045-7825(92)90042-I.
URL https://www.sciencedirect.com/science/article/pii/S004578259290042I

[43] T. Belytschko, M. O. Neal, Contact-impact by the pinball algorithm with penalty and lagrangian methods, International Journal for Numerical Methods in Engineering 31 (3) (1991) 547-572. arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1002/nme.1620310309, doi:https://doi.org/10.1002/nme.1620310309.
URL https://onlinelibrary.wiley.com/doi/abs/10.1002/nme.1620310309

[44] T. Belytschko, I. Yeh, The splitting pinball method for contact-impact problems, Computer Methods in Applied Mechanics and Engineering 105 (3) (1993) 375-393. doi:https://doi.org/10.1016/0045-7825(93)90064-5.
URL https://www.sciencedirect.com/science/article/pii/S0045782593900645

[45] M. Figueiredo, L. Marcelino, T. Fernando, A survey on collision detection techniques for virtual environments, Proc. of V Symposium in Virtual Reality, Brasil 307 (01 2002).

[46] C. Geuzaine, J.-F. Remacle, Gmsh: A 3-D finite element mesh generator with built-in pre- and post-processing facilities, International Journal for Numerical Methods in Engineering 79 (11) (2009) 1309-1331. doi:10.1002/nme.2579.
URL http://dx.doi.org/10.1002/nme.2579