Susanne Claus

Pierre Kerfriden

Faezeh Moshfeghifar

Sune Darkner

Kenny Erleben

Christian Wong

Contact modeling from images using cut finite element solvers

Keywords: Contact, Nonlinear Algorithm, LaTIn, nonconforming finite element, CutFEM, Ghost Penalty

teaching and research institutions in France or abroad, or from public or private research centers.

Introduction

In this article, we present a modern digital pipeline for the simulation of the mechanical interaction between solid bodies from CT images. The generation of finite element models from medical or engineering images is a very challenging task. The first step is to identify the set of voxels corresponding to the different solid bodies using image segmentation and/or registration. The regional information of the bodies is then typically transformed into a surface representation of each body, e.g. surface triangulations. From the simulation point of view, the challenge is that this surface mesh representation is of low quality and frequently contains badly shaped triangles which are not suitable for finite element analysis. Therefore, the surface mesh needs to be re-meshed, and this boundary operation is followed by interior volume meshing. Unfortunately, meshing is error-prone and often requires time consuming manual interventions and corrections, which is highly cumbersome, and requires valuable experts' time. To address this issue, several approaches are followed in the literature. A first class of mesh mapping methods [START_REF] Viceconti | A comparative study on different methods of automatic mesh generation of human femurs[END_REF] attempt to automatise mesh generation by deforming reference meshes to fit to specific geometries. Of course, such approaches still require prior meshing work and are limited 1 to relatively small variations around the reference geometry. A second class of approaches aim to eliminate or significantly reduce the reliance on meshing. The most popular instance of these methods is the voxel or image-based finite element method [START_REF] Keyak | Automated three-dimensional finite element modelling of bone: A new method[END_REF][START_REF] Luca Giovannelli | Direct medical image-based Finite Element modelling for patient-specific simulation of future implants[END_REF][START_REF] Ruess | The finite cell method for bone simulations: Verification and validation[END_REF]. In voxel based finite element methods, voxel of the image are directly converted to a hexahedral finite element mesh. Grey-scale values of the voxels are then mapped onto material values such as elasticity moduli. These methods reduce the meshing burden considerably, because the computational grid does not need to be aligned with any object boundary. However body boundaries are represented by "stair-cases". This causes stress singularities in the finite element simulation, which negatively affects the quality of the results near boundaries. Unfortunately, boundaries are often of key importance to analysts. This is the case in the present contribution, which focusses on the simulation of contact between two elastic bodies. In addition to this drawback of voxel-based methods, one may also find that, depending on the resolution of the image, two bodies may be indistinguishable and merged together into one big block. In this article, we introduce an alternative approach related to X-FEM [START_REF] Moës | Extended finite element method for cohesive crack growth[END_REF][START_REF] Fries | The intrinsic XFEM: A method for arbitrary discontinuities without additional unknowns[END_REF][START_REF] Legrain | High order X-FEM and levelsets for complex microstructures: Uncoupling geometry and approximation[END_REF] and CutFEM [START_REF] Burman | CutFEM: Discretizing geometry and partial differential equations[END_REF][START_REF] Hansbo | A finite element method for the simulation of strong and weak discontinuities in solid mechanics[END_REF]. We start from surface triangulations generated by a segmentation tool. This surface mesh, which may be of poor quality, is then embedded in a regular and fixed background grid. We now generate signed distance functions for each body, using the regular background grid to represent them discretely. As per usual with level-set-based methods, the surfaces of the bodies will be represented by the zero contour line of these functions. Zero contour lines may intersect the background mesh in an arbitrary manner. This allows for a sharp and smooth (continuous and elements-wise linear in this paper) surface representation which are then used simulate contact and compute contact forces between the elastic bodies. The background mesh which is used to discretize the unknowns and the level set functions needs to be sufficiently fine in order for the piecewise linear level set function to capture the geometry accurately. For this purpose, we develop and analyse two refinement criteria for adaptive mesh refinement. The contact problem between elastic bodies will be solved using a novel CutFEM-LaTIn algorithm first introduced in [START_REF] Claus | A stable and optimally convergent LaTIn-CutFEM algorithm for multiple unilateral contact problems[END_REF], but applied so far to mathematically defined geometries. We previously illustrated the versatility of this algorithm, and showed that it converges optimally. As an additional novel contribution in this article, we formally show that our contact algorithm is an extension of the augmented Lagrangian CutFEM algorithm [START_REF] Burman | Deriving robust unfitted finite element methods from augmented Lagrangian formulations[END_REF][START_REF] Burman | Augmented Lagrangian finite element methods for contact problems[END_REF]. These three elements, an efficient surface triangulation to level-set geometrical algorithm, adaptive refinement criteria and a state-of-the art contact simulation algorithm for contact over implicitly defined geometries, give rise to a stable and robust digital pipeline from CT images to finite element simulations.

This article is organised as follows. In Section 2, we detail the CutFEM-LaTIn algorithm used to compute the contact forces between two elastic bodies. In Section 3, we describe our algorithm to convert CT-scans to cut finite element biomechanical simulations for multiple bones including two adaptive refinement criteria. In Section 4, we show numerical results for contact between a femur and a hip bone.

Frictionless unilateral contact between two deformable solids

In this Section, we will describe our contact algorithm between two elastic bodies which occupy the domain Ω 1 and Ω 2 , respectively. The domains will be described implicitly using two level set functions. We begin our problem definition in the continuous setting in this subsection, before we describe its discretization of the domains and the definition of the discretization of the unknowns in the next subsection.

We assume that each elastic body is only undergoing small deformations and that they interact through frictionless unilateral contact. A contact problem between more than two bodies is possible in our framework and a contact algorithm for a larger number of solids has been introduced in [START_REF] Claus | A stable and optimally convergent LaTIn-CutFEM algorithm for multiple unilateral contact problems[END_REF]. However, in this contribution, we restrict ourselves to the contact problem between two bodies (see Figure ??). Let Ω b denote our total domain which contains our two elastic bodies. Let φ 1 : Ω b → R and φ 2 : Ω b → R be the two level set functions that describe the geometry of body 1 and the geometry of body 2, respectively. We then define two domains Ω 1 and Ω 2 by

Ω 1 := {x ∈ Ω b | φ 1 (x) ≤ 0 and φ 2 (x) > 0}, Ω 2 := {x ∈ Ω b | φ 2 (x) ≤ 0}. ( 1
)
The contact interface Γ c between the two solids, which is assumed to be of non-zero measure, is defined as

Γ c = {x ∈ Ω b | φ 2 (x) = 0 and φ 1 (x) ≤ 0}. (2)
Furthermore, let ∂Ω i denote the boundary of Ω i , i = 1, 2 and ∂Ω := (∂Ω 1 ∪ ∂Ω 2 )\Γ c the boundary of domain Ω 1 and Ω 2 without Γ c . We further partition the boundary ∂Ω into a Neumann part ∂Ω t and a Dirichlet part ∂Ω u with ∂Ω t ∩ ∂Ω u = ∅.

The problem of unilateral contact can be formulated as an optimisation problem under inequality constraints. We look for displacement field u

:= {u 1 , u 2 } : Ω 1 × Ω 2 → R 3 × R 3 that minimises the potential energy J (u) := 1 2 a(u, u) -l(u) (3)
under inequality constraints

u • n Γ ≤ 0 on Γ c (4)
and Dirichlet boundary conditions

u = u d on ∂Ω u . (5)
Here, n Γ denotes the normal on Γ c pointing from Ω 1 to Ω 2 , u d is a known prescribed displacement and u = u 1 -u 2 denotes the jump of u across the contact interface Γ c . The inequality constraint (4) enforces a zero interpenetration between the domains Ω 1 and Ω 2 . In the previous expression, the bilinear form a is defined as

a(u, δu) = a 1 (u 1 , δu 1 ) + a 2 (u 2 , δu 2 ) = 2 i=1 Ω i σ(u i ) : (δu i ) dΩ (6)
where σ(u i ) = λ i tr( (u i )) I I I + 2 µ i (u i ) is the Cauchy-stress tensor, (u i ) = 1 2 ∇u i + (∇u i ) T is the strain tensor, λ i and µ i are the two Lamé coefficients, i.e.

λ i = E i ν (1+ν)(1-2ν) , µ i = E i 2(1+ν)
, with Young's modulus E i and Poisson's ratio ν = 0.3, i = 1, 2. The linear form l is defined as

l(δu) = l 1 (δu 1 ) + l 2 (δu 2 ) = 2 i=1 Ω i f • δu i dΩ + 2 i=1 ∂Ω t ∩∂Ω i t d • δu i dΓ (7)
Here, f is a given body force and t d is a given boundary traction.

The primal constrained optimisation problem (3)-( 5) can be reformulated as an unconstrained primal/dual problem, which will be the basis for the development of this paper. In this setting, the solution to the contact problem is found by extremising the Lagrangian

L(u, λ) := J (u) + b(λ, u) (8)
with respect to primal variable u, in the space of fields satisfying the Dirichlet boundary conditions, but not the unilateral contact condition, which is now relaxed and enforced by

b(λ, u) = Γ c λ ( u • n Γ ) dΓ. (9)
We reformulate the primal dual formulation in terms of the Karush-Kuhn-Tucker conditions, which consist of the linear variational form

a(u, δu) + b(λ, δu) = l(δu) (10)
for all δu satisfying δu = 0 over ∂Ω u and the nonlinear constraint on the contact interface Γ c [START_REF] Dapogny | Computation of the signed distance function to a discrete contour on adapted triangulation[END_REF] which can be reformulated as

u • n Γ ≤ 0 λ ≥ 0, (u • n Γ ) λ = 0,
Γ c λ -γ u • n Γ + 1 γ λ + δλ dΓ = 0 (12)
for all δλ ∈ L 2 (Γ c ). Here, γ is a strictly positive scalar and [x] + denotes the positive part of a scalar quantity x ∈ R, defined as

[x] + = x if x > 0, 0 otherwise. ( 13
)
For a proof of the equivalence of condition ( 11) and (12) see for example [START_REF] Claus | CutFEM Method for Stefan-Signorini Problems with Application in Pulsed Laser Ablation[END_REF][START_REF] Chouly | A Nitsche-based method for unilateral contact problems: Numerical analysis[END_REF].

CutFEM primal/dual solver for the contact problem

In this Section, we will describe the discretization of the contact problem, which uses level set descriptions of the geometries in a fixed and regular background grid. The boundary of the elastic bodies will be represented by the zero contour lines of level set functions, which will be allowed to intersect the background mesh cells in an arbitrary manner. Each domain Ω i will be embedded in their own fictitious domain resulting in overlapping fictitious domains along the contact interface. Each u i will obtain its own set of degrees of freedom in its fictitious domain allowing for a multiple valued displacement u = (u 1 , u 2 ) along the interface. This enables jumps in the displacement along the interface.

CutFEM/Level-set approximation of the geometry

Let us introduce a triangulation T h of the background domain Ω b . We assume a simple shape for Ω b , which is easy to mesh with a regular background grid, e.g. a bounding box which contains both elastic bodies. Furthermore, let us introduce the finite element space of continuous piecewise linear functions, i.e.

Q h := {φ ∈ C 0 (Ω b ) : φ| K ∈ P 1 (K) ∀K ∈ T h }. (14)
We now define a piecewise linear approximation φ 1 h ∈ Q h of level set φ 1 and φ 2 h ∈ Q h of level set φ 2 , such that for any i ∈ {1, 2}, the nodal values of φ i h coincides with the values of φ i at the locations of the nodes. We now define the approximate physical domains Ω 1 h and Ω 2 h as follows:

Ω 1 h = {x ∈ Ω b | φ 1 h (x) ≤ 0 and φ 2 h (x) ≥ 0} Ω 2 h = {x ∈ Ω b | φ 2 h (x) ≤ 0} (15)
and the piecewise linear approximation of the interface

Γ c h = {x ∈ Ω b | φ 1 h (x) = 0 and φ 2 h (x) ≤ 0}. ( 16
)

Overlapping domain decomposition

Now, we define the set of all elements of T h that have a non-zero intersection with

Ω i h , i = 1, 2 T i h := {K ∈ T h : K ∩ Ω i h = ∅}. (17)
The domain corresponding to this set, denoted by Ωi h := K∈ T i h K, is called fictitious domain. Furthermore, let us denote all elements which are intersected with Γ c h by Ĝh :

= {K ∈ T h | K ∩ Γ h = ∅} (18)
The domain corresponding to this set is denoted by Γh := K∈ Ĝh K. This definition will be used to extend interface quantities to a band of elements.

Of key importance to CutFEM methods is the careful tailoring of regularization terms to ensure that the system matrix remains well-conditioned. To this end, we define the set of ghost-penalty element edges for fictitious domain Ωi

h , i ∈ {1, 2} as Fi G := {F = K ∩ K : K ∈ T i h and K ∈ T i h , F ∩ Γh = ∅}. (19)
For stabilization purposes still, we further define the set of intersected interface edges

FI := {F = K ∩ K : K ∈ Ĝh and K ∈ Ĝh }. (20)
Note our consistent use of notation . to indicate extended quantities.

Fictitious domain and extended interface finite element spaces

We seek the displacement approximation ûh = û1 h , û2 h of the two-body contact problem in the product space

Ûh = Û1 h × Û2 h , Ûi h := {û ∈ C 0 ( Ωi h ) : û| K ∈ P 1 (K) ∀K ∈ T i h }, i = 1, 2. ( 21
)
Note that ûh is multi-valued. This feature allows for the representation of embedded discontinuities in CutFEM. In all elements that are intersected by the contact interface Γ c h , û1 h and û2 h may take different values at a single point, which results in a non-zero displacement jump between û2 h and û1 h . Next, we define the space of extended interface fields of polynomial order 1 as

Ŵh := { Ŵ ∈ C 0 ( Ĝh ) : Ŵ | K ∈ P 1 (K) ∀K ∈ Ĝh }. ( 22
)
This is another important aspect of the CutFEM approach. Fields on manifolds such as embedded interfaces are represented as the trace of a finite element field defined in R d . In this paper, this feature will be used to represent the Lagrange multipliers required to enforce contact conditions in a primal/dual setting.

Primal/dual CutFEM solver

We now describe our P1/P1 LaTIn/CutFEM primal/dual algorithms, first introduced in [START_REF] Claus | A stable and optimally convergent LaTIn-CutFEM algorithm for multiple unilateral contact problems[END_REF]. We seek the displacement field in the linear CutFEM space Ûh and the Lagrange multipliers, which correspond to the normal reaction forces from contact, in the spaces of linear interface fields Ŵh . The presentation formalism that we have chosen for this paper allows us to highlight its relationship with the work presented in [START_REF] Burman | Deriving robust unfitted finite element methods from augmented Lagrangian formulations[END_REF][START_REF] Burman | Augmented Lagrangian finite element methods for contact problems[END_REF] in the specific case of frictionless unilateral contact (our algorithm extends straightforwardly to more general interface conditions such as friction contact or elastic damageable cohesive interfaces [START_REF] Kerfriden | A three-scale domain decomposition method for the 3D analysis of debonding in laminates[END_REF]).

Regularization of the bulk displacement

Let us define the extended and stabilized bilinear form over the product of fictitious domain fields as

âh (û h , δû h ) = 2 i=1 âi h (û h , δû h ) + 2 i=1 ŝi h (û h , δu h ) ( 23
)
where

for i ∈ {1, 2} âi h (û h , δû h ) = Ω i h σ(u i h ) : (δu i h ) dΩ (24)
and the so-called ghost-penalty terms are defined by

ŝi h (û h , δû h ) = F ∈ F i G F α h E i σ(u i h ) • n F • σ(δu i h ) • n F dΓ. (25)
Here, x • n F denotes the normal jump of the quantity x over the face, F , defined as

x • n F = x| K n F -x| K n F ,
where n F denotes a unit normal to the facet F with fixed but arbitrary orientation and α > 0 is the ghost penalty parameter.

Regularization of the Lagrange multiplier field

The discretized conditions for optimality are the following: Find (û h , λh ) ∈ Ûh × Ŵh satisfying the Dirichlet boundary conditions such that ∀(δ ûh , δ λh ) ∈ Ûh × Ŵh , where δ ûh satisfies the homogeneous boundary conditions (and there is no positivity requirement for the Lagrange multipliers),

       âh (û h , δû h ) + b h ( λh , δû h ) + ♥ s ♥ ( λh , δû h • n Γ h ) = l(δû h ) δ λh , λh -γ ûh • n Γ h + 1 γ λh + Γ c h + s ♥ ( δλ h , λh ) = 0 ( 26
)
where s ♥ is a regularization term that penalizes the jump of Lagrange multipliers across interface element edges, as introduced in [START_REF] Burman | Augmented Lagrangian finite element methods for contact problems[END_REF]:

s ♥ (δ λh , λh ) = βh 3 F ∈ FI F ∇δ λh • n F • ∇ λh • n F dΓ. (27)
Here, (u, v) Γ c h = Γ c h uv dΓ. And ♥ is either 0 or 1. The toggle boolean ♥ is introduced for comparison to other references. When ♥ = 1, the previous problem is equivalent to the extremization of the saddle-type functional

M(û h , λh ) := 1 2 âh (û h , ûh ) + γ 2 ûh • n Γ h + 1 γ λh + , ûh • n Γ h + 1 γ λh + Γ c h - 1 2γ b h λh , λh -l(û h ) -s ♥ ( λh , λh ) (28)
with respect to both of its arguments. However, in this contribution, we choose ♥ = 0.

Domain decomposition formulation

Whilst the previous "jump" formulation (26) couples the CutFEM subdomains, the LaTIn algorithm introduced in [START_REF] Claus | A stable and optimally convergent LaTIn-CutFEM algorithm for multiple unilateral contact problems[END_REF] implements a nonlinear domain decomposition approach. In this LaTIn-CutFEM algorithm the bulk problem is first solved in each subdomain independently and then interface fluxes and displacements are computed at the contact interface Γ c h . We iterate between the bulk problems (linear stage) for i = 1, 2 and the interface problem (local stage) until we achieve convergence. In the following, we briefly outline the main steps of the algorithm, a more detailed description can be found in [START_REF] Claus | A stable and optimally convergent LaTIn-CutFEM algorithm for multiple unilateral contact problems[END_REF]. To reformulate the "jump" formulation (26) into a problem over each fictitious domain Ωi h , we realize that ( 26) can be decomposed into

           2 i=1 âi h (û i h , δû i h ) + b h ( λh , δû i h ) = 2 i=1 l i (δû i h ), δ λh , λh -γ ûh • n Γ h + 1 γ λh + Γ c h + s ♥ ( δλ h , λh ) = 0, . (29)
where b h decomposes into b 1 h and b 2 h as

b h ( λh , δû h ) = Γ c h λh ( δû h • n Γ h ) dΓ = Γ c h λh (δû 1 h • n Γ h ) dΓ - Γ c h λh (δû 2 h • n Γ h ) dΓ =: b 1 h ( λh , δû 1 h ) -b 2 h ( λh , δû 2 h ) . (30)
Let k denote the iteration steps starting with k = 0. We then iterate between

The linear stage In which we solve the bulk problems: Find the displacement ûh ∈ Ûh such that for all i = 1, 2 and for all δ ûh ∈ Ûh

âi h (û i,k+1 h , δû i h ) + b i h ( λi,k+1 h , δû i h ) = l i (δû i h ) (31)
subject to the descent direction for the interface fluxes λi.k+1 h ∈ Ŵh and the interface displacements ŵi,k+1 ∈ Ŵh , i = 1, 2, ( λi,k+1

h - λi,k+ 1 2 h ) -γ( ŵi,k+1 -ŵi,k+ 1 2 ) = 0 , (32)
where the quantities from the previous half-iterate are known from the local stage and the interface displacements are given by ŵi,k+1 , δ λh

Γ c h = ûi,k+1 • n Γ h , δ λh Γ c h . ( 33
)
The local stage In each quadrature point along the interface Γ c h , find the interface fields { ŵi,k+

1 2 h , λ1,k+ 1 2 h }, i = 1, 2, that satisfy the contact condition λi,k+ 1 2 h = γ ( ŵ1,k+ 1 2 h - ŵ2,k+ 1 2 h ) • n Γ h + 1 γ λ1,k+ 1 2 h + (34) subject to search directions ( λi,k+ 1 2 h -λi,k ) + γ( ŵi,k+ 1 2 -ŵi,k ) = 0 , (35)
for all i = 1, 2 and to equality/definition λ2,k+

1 2 h = - λ1,k+ 1 2 h . ( 36
)
We regularize the interface quantities { ŵi,k+ , δ λh

Γ c h + s ♥ ( λi,k+ 1 2 h , δ λh ) = λi,k+ 1 2 h , δ λh Γ c h , (37)
for i = 1 or i = 2. The filtered interface jumps ŵi,k+ 1 2 for i = 1 and i = 2 can then be recovered locally using search direction (35), replacing the "double hat" -unfiltered -Lagrange multiplier by its regularized counterpart.

Digitial pipeline from CT-scans to CutFEM biomechanical simulations

In the following section, we will describe our pipeline from CT-images to CutFEM simulations. We will explain the pipeline on scans of the hip joint.

Image segmentation and surface triangulation of bone geometries

In the first step, we reconstruct a patient-specific hip joint geometry based on computed tomography (CT) images. The images are obtained from the cancer imaging archive (TCIA) [START_REF] Kirk | Radiology Data from The Cancer Genome Atlas Urothelial Bladder Carcinoma[END_REF] and transferred to the 3D slicer software for the hip joint segmentation [1,[START_REF] Fedorov | 3D Slicer as an Image Computing Platform for the Quantitative Imaging Network[END_REF]. We use multiple landmarks in a few slices to define the femur, the pelvis, and the background regions. Using these landmarks, all the image pixels are labelled and the corresponding 3D volumes are calculated (Figure 2). Further, we refine the bone models to eliminate rough surfaces, holes, and irrelevant connected tissues. We extract the surface of each bone as a triangulated surface mesh and export them as STL-files. These STL-files are used as the basis for the signed distance function calculation described in the following.

Proposed surface triangulation to level-set algorithm

The surface triangulations, read from the STL-files, are first embedded in a regular background grid and then signed distance functions are computed on this background grid. The following steps describe this surface triangulation to level set function algorithm.

Step 1: Determine Bounding Box Mesh In the first step, we determine a bounding box domain, denoted by Ω b , which contains both the surface triangulation of the hip bone, denoted by S 1 h , and the surface triangulation of the femur bone, denoted by S 2 h . This bounding box domain, Ω b , is defined by the minimum and maximum coordinates, x x x min = (x min , y min , z min ) and x x x max = (x max , y max , z max ), of all vertex coordinates in both surface triangulations S 1 h and S Step 2: Collision detection between the surface triangulations and the background mesh In the second step, we detect all collisions between the surface triangles and the background mesh tetrahedra. A collision between a surface triangle and a background mesh tetrahedron (cell) means that the surface triangle and the tetrahedron have a non-empty intersection. An efficient way to detect the collision between each surface triangle and the background mesh cells is to use Axis Aligned Bounding Box (AABB) trees. In this method, instead of testing each background mesh cell for an intersection with each surface triangle F j ∈ S i h , we detect intersections between two hierarchical structures of bounding boxes (AABB trees, for more details see [START_REF] Massing | Efficient implementation of finite element methods on nonmatching and overlapping meshes in three dimensions[END_REF][START_REF] Johansson | Multimesh finite element methods: Solving PDEs on multiple intersecting meshes[END_REF][START_REF] Ericson | Real-Time Collision Detection[END_REF]). This collision detection between an AABB tree for the surface mesh, S i h , i = 1, 2, and the AABB tree for the background mesh, T b h yields a set of tetrahedra that have a bounding box which collides with the bounding box around the surface triangle F j ∈ S i h . We then test among this set of tetrahedra with colliding bounding boxes if the surface triangle indeed intersects with the tetrahedra using geometric predicates (triangle -tetrahedra intersection [START_REF] Ericson | Real-Time Collision Detection[END_REF]). We express the results of this collision detection in a map between the index of each surface triangle, j for surface triangle F j ∈ S i h , i = 1, 2, and the indices of the background mesh tetrahedra, k m , it intersects (collides) with. Formally, these maps, one for each surface triangulation, are defined as

χ i : S i h → T b h , (38)
j → {k 1 , k 2 , . . . , k Nj }, (39)
where i = 1, 2, j denotes the index of the surface triangle F j ∈ S i h and runs from 1 to the number of surface triangles in S i h and k m , m = 1, . . . , N j , denotes the index of the background mesh tetrahedra that the surface triangle j intersects with. Note that, one surface triangle in general intersects several tetrahedra in the background mesh. The number of tetrahedra intersected by surface triangle with index j is denoted by N j . In addition to these mappings, we also construct the inverse map from the index of the background mesh tetrahedra to the index of the surface triangles that intersect the tetrahedra defined as

χ i,-1 : T b h → S i h , k → {j 1 , . . . , j N k }. (40)
Here, k denotes the index of the background mesh tetrahedron, and j m , m = 1, . . . , N k denotes the index of the surface mesh triangle F jm ∈ S i h that has a non-empty intersection with the background tetrahedron k and N k is the number of surface triangles which intersect the background mesh tetrahedron with index k. These two maps are crucial for the next steps.

Step 3: Decompose cells in the background mesh into interior cells, intersected cells and exterior cells We give each background mesh tetrahedron, a set of two marker values, M = (M 1 , M 2 ), M 1 ∈ {0, 1, 2}, M 2 ∈ {0, 1, 2} depending on if the tetrahedron is inside,intersected or outside by femur or hip bone surface triangulations, S i h , i = 1, 2. We start with all cells marked as M = (0, 0) and then mark the tetrahedra intersected by S i h using the inverse mapping χ i,-1 . From the mapping, we determine the indices of all background mesh tetrahedra that have a non-zero intersection with the surface triangulation S i h , i = 1, 2 and mark these cells with 1 (M 1 = 1 if tetrahedron intersects with cells in S 1 h and M 2 = 1 if tetrahedron intersects with cells in S 2 h ). The last step in the mesh marking is to determine all cells outside of the femur bone surface triangulation (M 1 = 2) and outside of the hip bone surface triangulation (M 2 = 2). To achieve this, we assume that the background mesh cells are clearly separated by a band of intersected background mesh elements into mesh cells inside the femur or hip bone and into mesh cells outside. Starting with a tetrahedron in a corner of the background mesh, for which we are sure that the cell is outside of both the hip and femur bone geometries, we mark this corner cell with M = (2, 2). We then find all cells connected to this starting tetrahedron and see if any of the connected tetrahedron is marked with M i = 1. If the neighbouring cell is marked with M i = 1, we leave the mark at M i = 1 if it is marked with M i = 0, we set M i = 2 and find all cells connected to the new cell we just marked with M i = 2. This way, we find all cells that are connected to the outside corner cell without crossing the band of elements which are marked as intersected. For more details on a mesh marking method for connected domains see [START_REF] Dapogny | Computation of the signed distance function to a discrete contour on adapted triangulation[END_REF]. Note that the presented technique is similar to the classical coloring technique employed in Delaunay mesh generation [START_REF] Cheng | Delaunay Mesh Generation[END_REF]. Alternatively to exploiting connectivities of inside and outside domains, ray casting techniques can be used in which the number of intersections of a ray originating in a background mesh cell with the surface mesh can determine if a cell is inside or outside (see e.g. [START_REF] Baumgärtner | A robust algorithm for implicit description of immersed geometries within a background mesh[END_REF]). The decomposition into intersected, interior and exterior cells will be crucial in the next steps.

Step 4: Extract ficititious domain In this step, we disregard the cells outside the femur and hip bone geometries. The cells outside of both the femur and the hip bone have no physical meaning and are not of interest in our simulations. We disregard the outside cells by constructing a new background mesh, T * h ⊂ T b h , with all cells from the original background mesh which are marked with M 1 = 0 or M 2 = 0 and M 1 = 1 or M 2 = 1. In extracting the cells, we keep track of the mapping of the cell indices of this new "submesh" to the original background mesh to update the mappings χ i,-1 (to map from the surface triangle index to the new indices of the fictitious domain background mesh, T * h ). We also keep the mesh cell marker. The new background mesh, T * h , forms a "Cube"-like mesh which fully contains the hip joint geometry (as shown in Figure 3).

Step 2

Step 1

Step 4 Step 3 Step 5: Compute the signed distance functions to the femur and to the hip bone surface meshes. In this step, we determine two signed distance functions: one to the hip surface mesh, denoted by φ 1 h , and one to the femur bone surface mesh denoted by φ 2 h . We compute the distance to surface mesh S i h in the following way 1. Initialization phase. In the first step, we use a geometry based approach to compute the unsigned distance of all nodes of the intersected tetrahedra to the surface mesh, as described in the following. Using the mapping χ i,-1 , we iterate over all cells that are intersected by the surface mesh S i h . For each intersected cell, we iterate over its vertex nodes with coordinates v j p , j = 1, 2, 3, 4. We compute the distance from each vertex node to all surface triangles, which the current element intersects (given by map χ i,-1 ). This is done by calculating the unsigned distance of a point (v j p ) to a triangle [START_REF] Ericson | Real-Time Collision Detection[END_REF]. Among all distances determined for a vertex node, we keep the shortest distance and save the value in the degree of freedom of the level set finite element function corresponding to the vertex node given by a vertex to dof mapping (for a piecewise linear level set function). Note here, that vertex nodes are often iterated over several times (when they belong to several intersected tetrahedra) and we always keep the shortest unsigned distance value of all computed distances.

2. After initialising all values of the level set at the nodes of all intersected tetrahedra (previous step), we use a fast marching method [START_REF] Sethian | Fast marching methods[END_REF] to compute the approximate unsigned distance of all remaining nodes in the background mesh. In a first step, we iterate over all nodes with known distance value and determine all nodes connected to these nodes via an edge (neighbouring node). We save these connected (neighbouring) nodes in a set called active set. We then compute the distance of these nodes in the active set to the known values, as follows. For each node in the active set, we iterate over all tetrahedra that share this node. For each tetrahedron, we compute the distance of the active node to known values in the tetrahedron (if any). There are four possible cases, firstly, the tetrahedron contains three known values (point (active node)triangle distance), secondly, the tetrahedron contains two known values (point (active node)edge (between known nodes) distance); thirdly, the tetrahedron contains one known value (point (active node) -point (known node) distance) or the tetrahedron contains no known values (no distance computation). Once all nodes in the active set have their distance values (smallest of all computed distances), we determine the next set of nodes connected to all known nodes. We repeat this process until all unsigned distance values in all nodes have been computed. For more details on the employed fast marching method see [START_REF] Reusken | Numerical Methods for Two-Phase Incompressible Flows[END_REF].

3. In a last step, we convert the unsigned distance function computed in the previous steps to a signed distance function using the cell marker computed in step 3. In step 3, we have marked all tetrahedra enclosed by the surface mesh with 0, all cells outside with 2 and all intersected cells with 1. We want negative level set values inside our bone geometries and a positive sign outside such that the zero contour of the level set function represents a linear approximation of the surface mesh (see Figure 4) and the normalised gradient of the level set represents a normal pointing from inside the bone to the outside. To achieve this, all that is needed is to assign negative signs to all values of the level set in nodes for cells marked with 0.

Following steps 5.1 to 5.3, for each S i h , i = 1, 2, we obtain the piecewise linear signed distance functions φ 1 h and φ 2 h , which provide a linear approximation of the femur and hip bone geometries, respectively. In the next Section, we will use these signed distance functions to compute stress distributions in the bones.

Adaptation of the background mesh to guarantee sufficient geometric resolution

As described in the previous section, we compute a piecewise linear signed distance function to represent the surface triangulation geometry. This piecewise linear representation is close to the surface geometry in planar areas. However, in areas of high curvature or in areas with fine features, we need to ensure that the background mesh size is fine enough to capture these features with our piecewise linear approximation. In this section, we will first make some heuristic observations about the quality of the level set approximation for the hip and femur bone geometry for a refinement of all background mesh cells which are colliding with the surface triangulations of the two bones. Secondly, we will use these heuristic observations to suggest two refinement criteria to automatically choose only some of the colliding elements to be refined and to generate a level set approximation that captures the geometry with sufficient accuracy. Firstly, Figure 5 (a) and (b) show the effect of choosing a background mesh, which is too coarse to represent the fine features in the hip bone geometry. The piecewise linear level set approximation cannot capture the fine bone structures in the upper hip bone and in the acetabulum and instead these areas are represented as holes in the geometry. Refining the mesh in the elements, which are colliding with both surface triangulation (as identified by the markers in the algorithm introduced in the previous section) until all features are sufficiently represented, can solve this problem. Figure 5, shows the zero contour line of the level set function for two such background mesh refinement steps (a,b,c) together with their background mesh (e,f,g) and the original surface triangulation (d,h). We can see in Figure 5 that the level set approximation leaves large holes in the hip bone for the initial background mesh. After the first refinement step, we still observe some small holes in the thin part of the upper hip bone which is finally completely closed in the second refinement step. The appearance of these holes in the hip bone geometry are due to sections of the hip bone geometry with no connectivity to interior cells as displayed in Figure 6. Figure 6 (a,b,c) shows a slice (shown in (g)) through the hip bone through a part of the bone with the holes. Elements which are colliding with the surface triangulation of the hip bone or the femur bone are displayed in green, elements which are outside the hip bone are shown in red and elements fully inside the hip bone are shown in blue. The zero level set contour line is shown in yellow and the surface triangulation is shown by spheres in the nodes of the surface triangulation. We observe that in the thin part of the hip bone, elements are detected as colliding but there is no element which is fully inside the hip bone geometry (see Figure 6 (d,e,f)). The computed distance function to the surface triangulation only gets assigned a negative sign for all nodes of elements that are fully inside. As the thin section does not contain any elements that are fully inside, the level set function does not change sign and hence there is no zero contour intersecting these elements. This yields the holes. By refining all background elements which are colliding with the surface triangulations, we can close these holes and improve the geometrical approximation. We refine until the zero level set contour is close to the original surface geometry.

Here, two refinements were found to be sufficient. As the femur bone does not contain so many small features, it is much easier approximated by a piecewise linear level set approximation in a coarse background mesh (see Figure 7). Figure 7 displays the linear approximation of the zero contour line of the level set function in yellow and the nodes of the surface triangulation as spheres in a slice through the femur bone. We observe that the zero contour line approximates the femur shape in the slice very well. Only for the coarsest mesh, we can observe a clear difference between the femur bone and the surface triangulation in areas of high curvature. Hence, we heuristically observed two major challenges in the resolution of the surface triangulation: the appearance of holes and the approximation of areas with high curvature. We develop two refinement criteria to address these two challenges. One to identify cells, in which the level set approximation does not change its sign even though the surface triangulation collides with the element (e.g. holes) and one to identify cells in which the zero contour line of the level set is far away from the surface triangulation (e.g. high curvature).

The first refinement criterion is to mark all cells for refinement which are colliding with the surface triangulation but in which the computed level-set function does not change sign. We refer to this refinement criterion as "adaptive sign criterion".

For the second refinement criterion, we determine for each background mesh cell which collides with at least one surface triangulation what the maximum absolute level set value is at all nodes of all surface triangles which collide with the background cell. For example, let us assume background mesh cell j collides with surface triangle k,l and m. We then evaluate the level set function value in each node of the surface triangles k,l and m, which will give us the distance of each node of the surface triangles to the zero level set contour. We take the maximum absolute level-set/distance value of all these node values and store this value for cell j, denoted by d j max . We repeat this process for all cells that collide with surface triangles. We then determine the global maximum value among all these cells, denoted by d global max ,

d global max = max j d j max ( 41
)
where j is a list of cell indices of cells that collide with the surface triangulation. We mark all cells for refinement which have a maximum distance value larger than a threshold times the global maximum, with 0 < < 1, i.e.

d j max > d global max ( 42
)
for each cell j which collides with the surface triangulation. We refer to the second refinement criterion as "adaptive greedy criterion".

To evaluate the quality of the level set approximation with respect to the surface triangulation, we calculate the average distance value, d avg , over all nodes of the surface triangulation, i.e.

d avg = k |φ i (x k x k x k )| N k (43)
where x x x k denotes the coordinate of a node in the surface triangulation, k = [0, 1, ..., N k ] is the node index, and N k is the total number of nodes in the surface triangulation. Figure 8 shows the average distance value of the level set function to the surface triangulation versus the number of cells for the different adaptive refinement criteria. We compare the refinment of all cells that are colliding with the interface (dashed, red) with the sign refinement criterion (dashed, cyan), the greedy refinement criterion with a parameter of = 0.5 and the combination of both sign and greedy criteria. Figure 8 (a) shows the average distance of the level set function to the surface triangulation for the femur bone and Figure 8 (b) shows the average distance for the hip bone. For each bone the refinement criteria have very different effects. For the smooth femur bone, the sign criterion, which will refine in the area of holes or very high curvature does not refine as efficiently as greedy or uniform refinement with respect to the average distance. For the hip bone on the other hand, the sign criterion reduces the average distance between the level set function and the surface triangulation most efficiently. This is what we have heuristically observed for the hip bone before. The hip bone has a lot of fine feature and areas of high curvature and this is efficiently resolved by the sign criterion. The greedy criterion performs best for the femur bone and yields lower average distances than uniform refinement for the same number of elements. The greedy algorithm for the hip bone performs slightly less well for the hip bone than the sign criterion up until 70000 cells by which the greedy criterion yields better refinement than the sign criterion. The combination of both sign and greedy criterion yields the lowest average distance to the surface triangulation for the hip bone geometry.

Unilateral contact between hip and femur bone

In this section, we use our LaTIn/CutFEM algorithm described in section 2 to compute the stress distribution in the hip and the femur bone in unilateral contact. All of our software developments are based on FEniCS [START_REF] Alnaes | The fenics project version 1.5[END_REF]. In [START_REF] Claus | A stable and optimally convergent LaTIn-CutFEM algorithm for multiple unilateral contact problems[END_REF], we have shown that our LaTIn/CutFEM algorithm was stable and converged optimally for a range of carefully selected benchmark problems. This contribution focuses on showing that our contact algorithm performs well for geometries obtained from images. In this section, we use the level set geometry description for the femur and the hip bone obtained from the algorithm described in the previous section and compute a unilateral contact problem on the bones. We bring the femur and the hip bone into contact by moving the femur bone geometry slightly into the acetabulum to create a clear contact interface in between the bones. The resulting contact interface is displayed in Figure 9. Here, we need to predefine the contact interface as our algorithm currently does not include a contact detection algorithm and our algorithm is only valid for small deformations.

The extension of the algorithm to include large deformations and contact detection will be addressed in a future publication but is out of the scope of this contribution.

Figure 10 shows the fictitious domain mesh used for the computation, which consists of 460.985 cells yielding 288.426 number of degrees of freedom for the P1 displacement. As boundary conditions, we apply homogenous Dirichlet boundary conditions for the displacement on the bottom of the femur bone and we apply a displacement of u = (0, 0, -1) T on the boundary facets marked in red in Figure 10 on top of the hip bone. We set the Young's moduli for the hip and femur bone to E 1 = E 2 = 1, and the penalty parameters to α = β = 0.1, γ = 1. Figure 11 shows the stress component σ zz for the hip and femur bone. We observe that under our artificial boundary conditions, the thin part of the hip bone is under compression and the rest of the hip bone is under tension. The femur bone is mainly under tension except for a small area in the femoral neck and in the contact region. Looking at a slice through the hip joint as shown in Figure 12 contact forces transmit the compression force to the top of the femoral head. Figure 12 (b) shows the deformation from the slice in the reference configuration to the slice in the deformed configuration. We observe that the femur bone deforms very little while the hip bone shows a much larger deformation. The magnitude of the deformation is shown in Figure 13, like seen in the slice the femur bone deforms very little while the hip bone deforms significantly away from the area of contact between the femur and hip bone. The femur bone stabilises the hip bone in the contact region but the hip bone is free to move away from the contact interface.

Conclusions

In this article, we have first presented a reformulation of our LaTIn/CutFEM algorithm in terms of an augmented Lagrangian approach for two elastic bodies. We have then detailed our robust algorithm to compute level set representations on finite element meshes for multiple bodies from segmented images.

We have discussed in detail on how to adapt the finite element mesh to guarantee a good approximation of the level set function to the surface triangulation. We have developed and analysed two refinement criteria for adaptive mesh refinement. Finally, we have concluded our article with numerical results of a unilateral contact computation between a femur and a hip bone. Our contact algorithm is currently limited to small deformation. In a future contribution, we are planning to extend our approach to large deformation.

Figure 1 :

 1 Figure 1: Contact problem between two domains Ω 1 and Ω 2 described via two level set functions φ 1 and φ 2 and contact interface Γ c embedded in background domain Ω b .

  2 h . The resulting domain is given by Ω b = [x min , x max ] × [y min , y max ] × [z min , z max ]. We then generate a uniform and regular tetrahedral mesh, denoted by T b h , on the background box domain Ω b .

Figure 2 :

 2 Figure 2: The different scene views of the hip joint and the calculated 3D model of both the femur (green) and the pelvis (yellow) bones in 3D slicer [1].

Figure 3 :

 3 Figure 3: Illustration of the first 4 steps in the surface triangulation to level set algorithm. The surface triangulation is first embedded in a regular background grid (Step 2). The cells are then marked depending on if they are inside, intersected or outside of the femur and hip bone geometries (Step 3). All cells which are inside or intersected by either the femur or the hip bone geometry are extracted (Step 4).

Figure 4 :

 4 Figure 4: Linear approximation of the geometry (zero level set contour) given a surface triangulation for the femur bone. The grey balls indicate the nodes of the original surface triangulation. All elements with positive level set values in all their nodes are shown in red, all elements with a change in sign of the level set value (containing the zero contour line of the level set function) are shown in green and all elements with negative level set values in all their nodes are shown in blue.

Figure 5 :

 5 Figure 5: Effect of background mesh size onto geometric approximation. The zero level set contour for 2 steps of refinement of all tetrahedrons which are colliding with the surface triangulation is shown in (a),(b) and (c) together with their corresponding background mesh (e,f,g). The original geometry (d) together with the surface triangulation (h).

Figure 6 :

 6 Figure 6: Linear approximation of the hip bone surface triangulation in regular background mesh with two refinement steps of all elements which collide with the surface triangulations of the hip and femur bone. The figures (a), (b), (c) show elements marked as interior as blue, marked as colliding with green and marked as outside with red. The figures (d),(e),(f) shows all the elements which are fully inside the hip bone geometry.

Figure 7 :

 7 Figure7: Linear approximation of the femur bone surface triangulation in regular background mesh for two refinement steps in a slice through the femur head (a),(b), (c) and a zoom to the area of high curvature in the femur bone slice (d),(e),(f). For the coarsest mesh, we observe a clear difference between the piecewise linear approximation of the level set function and the surface triangulation in the area of high curvature. This difference is resolved by mesh refinement.

Figure 8 :

 8 Figure 8: Comparison different adaptive refinement strategies in terms of average distance value for the femur bone level set function (a) and the hip bone level set function (b).

Figure 9 :

 9 Figure 9: Contact interface.

Figure 10 :

 10 Figure 10: Mesh and boundary conditions.

Figure 11 :

 11 Figure 11: Stress component σ zz in the hip and femur bone in front view (a) and side view (b).

Figure 12 :

 12 Figure 12: Stress component σ zz (b) in slice of hip joint (a) together with its deformation (b).

20

 20

Figure 13 :

 13 Figure 13: Magnitude of deformation in hip joint.

Funding The authors gratefully acknowledge the financial support provided by the European Unions Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No. 764644 (ITN RAINBOW) and the funding received by the Alfred Benzon Fond. The publishing fees of this article have been graciously covered by the French Association for Computational Mechanics (CSMA).

Declarations

Availability of data and materials. The CT-scan of the hip joint used in this article was obtained from the cancer imaging archive (TCIA) [START_REF] Kirk | Radiology Data from The Cancer Genome Atlas Urothelial Bladder Carcinoma[END_REF] and can be downloaded at https://www.cancerimagingarchive.net/.

Competing interests. The authors declare that there are no competing interests.

Authors' contributions. All authors contributed equally to this article.

Acknowledgements. Not applicable.