1 Introduction
Unfitted or nonconforming finite element methods uncouple the description of the geometry from the representation of the solution field itself. Typically, the geometry of the computational domain is projected over a regular background grid. In this setting, boundaries or interfaces between objects cut through elements of the corresponding mesh. Nonconforming methods are attractive for applications where coupling analysis codes and third party meshing libraries is either impractical, numerically expensive and/or prone to errors. Nonetheless, a number of specific challenges needs to be addressed for unfitted numerical solutions to be computable. Firstly, integrals need to be calculated over cut elements, which requires specialised numerical quadratures. Secondly, unfitted approaches require stabilisation. This is because combinations of degrees of freedom may be poorly controlled in regions where contributions from cut elements to integral forms are small. Finally, in the context of multiple interacting materials, enrichment of the finite element solution space is required in order to allow for numerical jumps or kinks to develop over embedded interfaces. Failure to do so may severely impair the convergence rate of the unfitted finite element solver.
Over the last two decades, several encompassing frameworks have been developed to provide guidance for the development of unfitted finite element solvers. The eXtended Finite Element Method (XFEM) [36, 18, 7, 21] relies on the Partition of Unity Method [34] to enrich the approximation space. Within this popular framework, nonconforming discontinuities can be introduced in an elegant and robust manner. Field discontinuities can also be enabled by making use of Embedded Discontinuity Methods [39, 37, 31], whereby discontinuous shape functions are directly defined at the element level. The CutFEM framework [22, 10], which is fundamental to the developments reported in this paper, uses an alternative concept to introduce nonconforming field discontinuities, and to perform computations involving nonconforming geometrical domains. In the context of multiple interacting materials, CutFEM solvers enrich the solution space by overlaying (i.e. doubling or tripling) cut elements. Subsequently, penalty techniques such as the Nitsche method are employed to disallow unwanted discontinuities, for instance to represent jumps in fluxes whilst preserving the continuity of the primal field [22]. CutFEM solvers are usually equipped with a ghostpenalty stabilisation technique [9], which results in optimal convergence rates whilst ensuring that the condition number of the system matrix degrades at the same rate as that of the conforming finite element method when the mesh is refined.
Extending unfitted finite element methods to nonlinear problems, and especially to problems involving complex interface behaviour, remains a challenge that needs to be addressed on a casebycase basis. The present paper focuses on problems involving unilateral contact between elastic solids. This class of problems is important to engineers, and they are typically encountered when modelling joints between components in largescale engineering assemblies. Numerically, unilateral contact is notoriously difficulty to treat, even in the context of fitted finite element approaches. References on the matter include [45, 30, 42, 1, 33]. Of particular interest to us are LaTInbased solvers [13, 28, 8, 32, 29, 2], which treat unilateral contact through a twofield formulation of the interface equations, and resolve these fields iteratively using a dedicated twostage solver. A number of researchers have demonstrated the versatility of LaTIn solvers, owing to the fact that the predictor stage (socalled “linear stage”), is weakly dependent on the type of interface conditions, and completely independent of the current interface state. As a result, switching between nonlinear interface conditions is relatively painless [24]. Moreover, dedicated spacetime Model Order Reduction techniques may take advantage of the predictor invariance to reduce the cost of timedependent or multiquery computations [29, 24, 19]. However, a noteworthy issue related to the hybridmixed LaTIn formulation is the difficulty of discretising displacement and traction interface fields in a consistent manner. Whilst a piecewise constant discretisation of interface fields seems to be physically sound when associated with P1 finite element spaces for the bulk equations, such a discretisation scheme allows for the onset of numerically uncontrolled deformation modes whose associated amplitude creep up with increasing LaTIn iteration count, and eventually alter the convergence rate of the finite element solver. A claimed remedy to this numerical plague is to introduce additional degrees of freedom for the bulk fields, by means of either hrefinement or prefinement. An alternative to this numerically expensive procedure has recently been proposed in [20], as discussed later on. Although both approaches seem to regularise the contact problem efficiently, the effect on the convergence of the mixed finite element solver with mesh refinement has not been thoroughly studied.
Over the years, several attempts to develop stable unfitted finite element solvers for unilateral contact have been met with success. XFEMbased solutions have been proposed in [16, 25, 40, 17, 32, 20, 38]. Of particular importance to us are LaTInbased unfitted solvers, as the LaTIn is also a fundamental technological element for the present contribution. The first of these solvers was proposed in [16] to solve crack propagation problems (i.e. the historically favoured playground for XFEM). In [40], the XFEMLaTIn idea was extended in order to simulate fatigue crack propagation. In [32], the authors compared a LaTIn formulation to a penalty/Newton approach, and claimed the superiority of the latter in terms of convergence rate and stability. Finally, the authors of [20] proposed a regularised, LaTIntype formulation of contact for XFEM. The main idea of the approach is to relax the kinematic continuity condition between the bulk primal field and its interface trace and to enforce it weakly via a regularised mortar method.
In contrast to XFEMbased developments, the CutFEM framework is relatively underdeveloped in the context of unfitted unilateral contact. We mention the recent work described in [14, 15], which extends the Nitschetype treatment of interface conditions that is usually employed in CutFEM solvers to the case of unilateral contact. The resulting fitted finite element solver is proven to be optimally convergent, stable and numerically efficient. A very first nonconforming CutFEM strategy for unilateral contact is discussed in [11], whereby the contact conditions are treated by an augmented Lagrange multiplier formulation, and the discrete equations are stabilised by a ghostpenalty regularisation.
The present paper proposes a novel LaTInCutFEM solver for unilateral contact between multiple interacting solid bodies. The solid bodies will be described by linearised elasticity, which implies that the location of contact interfaces is known a priori. Our algorithm relies on the traditional CutFEM enrichment strategy [22]. However, contact conditions over nonconforming interfaces will be enforced via the LaTIn method. The resulting algorithm is a twostage solver, whereby the linear stage (i.e. linear predictor) consists in solving a series of independent linear problems for each of the solid bodies, whilst the local stage (i.e. interface corrector) is nonlinear. The method naturally inherits the coarsegrain parallel characteristics of LaTInbased domain decomposition methods [27, 28, 24], together with the previously mentioned versatility of LaTIn solvers with regards to the type of interface conditions being dealt with (i.e. unilateral contact with Coulomb friction [8, 12], frictional sliding [5, 6], cohesive fracture models with transition to unilateral frictionless contact [24]).
Two algorithmic elements of the proposed solver require particular care. Firstly, the LaTInCutFEM scheme needs to be stabilised. We regularise each of the subproblems of the linear stage using a ghost penalty technique [9]. We will show that, when using the appropriate scaling for the ghost penalty terms, the condition number of the system matrices “degrades optimally” with mesh refinement (i.e. it degrades at the same rate as that of the conforming finite element method). Secondly, we need to be careful when discretising the contact interface fields. As mentioned before, a naive discretisation of the LaTIn hybridmixed formulation leads to the onset of spurious interface solutions. We will show that the P1P1 discretisation scheme that we propose (i.e. P1 for the bulk primal field, P1 for both its trace and the trace of its normal flux) yields optimal convergence rates with mesh refinement, and eliminates the onset of instabilities. Our scheme “nonlocalises” the local stage of the LaTIn, as it relies on a twoscale treatment of the interface fields and nonlinear equations. We will also show that the nonlinear solver of the LaTIn converges fast, in the sense that it quickly leads to iterates whereby the discretisation error is larger than the algorithmic error, typically within twenty iterations. To our best knowledge, these convergence results have never been reported in the context of LaTIn solvers, even in the context of conforming geometries.
Our formal developments are accompanied by a highperformance computer implementation. The core of our implementation is the finite element C++/Python library FEniCS [4], which, in particular, proposes a range of highlevel tools to rapidly developed finite element solvers. The CutFEM C++ library, LibCutFEM, developed by S. Claus and A. Massing, and partially described in [10], defines additional tools that are specific to unfitted finite elements. This library forms the basis for the LaTInCutFEM code. We will present results in two and three space dimensions, which is made seamless by the FEniCS framework. Some of our examples exhibit a high level of geometric complexity, which is handled through a robust implementation of multiple level sets. In particular, we will study the convergence of contact problems exhibiting multiple interacting bodies (i.e. more than two) at a single point, (see other contributions on the topic, in a different context, in [5, 6]).
This paper is organised as follows. The LaTIn hybridmixed formulation of contact is presented in a continuous setting in Section 2, together with the associated twostage nonlinear algorithm. The nonconforming regularised LaTInCutFEM formulation is presented in Section 3. Section 4 is dedicated to numerical investigations, where we demonstrate the optimal convergence and stability properties of the proposed framework. We also showcase the versatility of our highperformance computer implementation.
2 LaTIn method for unilateral contact: the continuous setting
2.1 Domain decomposition
We consider a problem of linear elasticity over a 2 or 3dimensional spatial domain composed of nonoverlapping subdomains , where . The domain decomposition satisfies . Let us denote the set of all subdomains by . Subdomains may be multiplyconnected. The boundary of subdomain is denoted by , for any , whilst denotes the boundary of computational domain . Furthermore, we define interface domains for any such that the two indexes satisfy and the result of intersection is of nonzero measure. The union of these interfaces is denoted by . The set of pairs of indexes such that is defined as a valid interface is denoted by , and the corresponding set of interfaces is denotes by . We also denote the union of all the interfaces of subdomain by , for any subdomain index . We denote the outward normal pointing from to on by . We further note that we have , and we denote the outward normal to over by .
As an illustration for this set of definitions, consider the case of the multiphase domain depicted in Figure 1. The cubic domain comprises three subdomains , and . There are three interfaces between the domains: , and . Therefore, the interface of individual subdomains reads as , and .
As further detailed in the next subsection, subdomain field variables will be governed by the equations of linear elasticity, while contact equations to be solved over interfaces will couple subdomain problems to one another. The LaTIn framework [26] is particularly well suited to the treatment of such problems where complex nonlinear equations need to be solved over embedded interfaces [8, 29, 3]. In this particular context, the LaTIn method relies on two main ingredients:

a hybrid description of bulk and interface fields, and a mixed representation of the interface unknowns,

a twostep iterative solution algorithm for the solution of the global nonlinear problem formulated in its hybridmixed form.
Both these steps are described below.
2.2 HybridMixed LaTIn formulation of unilateral contact problems
Hybrid fields
For every subdomain , where , we introduce a displacement field , where (2D problem) or (3D problem).
For every domain interface with , we define two field variables:

a force field that represents a density of tractions (i.e. flux) applied to subdomain by adjacent subdomains through interface .

a displacement field that represents the trace of over .
As detailed below, will be governed by linear elasticity, while over an interface between two domains and , where , interface models will provide a link between the restriction of fields , to and the restriction of fields , to .
Bulk equations and boundary data
For every subdomain , , the displacement field is required to be sufficiently regular (continuous in particular), and to satisfy the following equations:

Static equilibrium equation and constitutive relation:
For all , the displacement field satisfies(1) In equation (1), is a body force and is the outer normal of domain . The stress function is defined as , where the strain is the symmetric part of the displacement gradient, i.e.