An insertability constraint for shape optimization

Patient-specific implants offer a host of benefits over their generic counterparts. Nonetheless, the design and optimization of these components present several technical challenges, among them being the need to ensure their insertability into the host bone tissue. This presents a significant challenge due to the tight-fitting nature of the bone-implant interface. This paper presents a novel insertability metric designed to efficiently assess whether a rigid body can be inserted into a tight-fitting cavity, without interference. In contrast to existing solutions, the metric is fully differentiable and can be incorporated as a design constraint into shape optimization routines. By exploiting the tight-fitting condition, the problem of planning an interference-free insertion path is reformulated as the search for a single interference-free movement, starting from the inserted configuration. We prove that if there exists any outward movement for which no interference is indicated, then the body can be fully extracted from or, equivalently, inserted into the cavity. This formulation is extremely efficient and highly robust with respect to the complexity of the geometry. We demonstrate the effectiveness and efficiency of our method by applying it to the optimization of two-dimensional (2D) and three-dimensional (3D) designs for insertability, subject to various design requirements. We then incorporate the proposed metric into the optimization of an acetabular cup used in total hip replacement (THR) surgery where geometric and structural requirements are considered.


Introduction
Insertability analysis, which is the problem of assessing whether a rigid body can be inserted into a tight-fitting cavity, is ubiquitous in a wide range of engineering applications.In automated tooling, for example, it is used to plan an interference-free path for inserting robotically guided parts into mating fixtures (Canny 1988;Oliver and Huang 1994).In molding, it is used to ensure that the part can be extracted from the core and cavity without excessive force (Carley 1993).And in engineering design, it is used to design tightly fitting parts and establish tolerances for assemblability (Shen et al. 2005).In particular, with respect to the design of orthopedic implants, where a tight fit with the host bone is usually required, insertability analysis plays a key role in avoiding obstructing geometry and stuck configurations.
With advances in medical imaging and computational analysis, virtual models of a patient's specific anatomy can be synthesized and used to develop custom orthopedic implants, designed to best fit the patient's unique anatomy.These implants may also be optimized for maximum longterm performance by simulating their effect on the host bony tissue.Nonetheless, the design of patient-specific implants is fraught with new technical challenges, not least of which is the requirement that the implant be insertable into the host bone tissue, without interference.The insertability requirement is particularly difficult to incorporate into a design optimization framework because of (1) the tight-fitting nature of the bone-implant system, and (2) the need for a differentiable insertability metric.
The problem of assessing the insertability of tight-fitting components has been extensively studied in the field Responsible Editor: Helder C. Rodrigues of robotics, where it is known as the peg-in-hole problem (Sacks and Joskowicz 1995).It is typically treated as a special instance of the general path planning problem, wherein the goal is to determine how to insert a movable body (the peg) into a stationary cavity (the hole) without interference (Latombe 2012).What makes the peg-in-hole problem challenging is that the body and cavity share nearly complementary geometry in the inserted configuration, resulting in highly constrained movement.This feature significantly decreases the efficiency and efficacy of traditional path-planning strategies, especially when considering 3D motion and complex geometry.Traditionally, solving the peg-in-hole problem involves searching the space of all possible physical configurations (configuration space) for a continuous non-overlapping set of interference-free configurations leading from an initial uninserted configuration to a final inserted configuration (Joskowicz and Taylor 1996).The search for an insertion path involves two tasks: (1) determining if a particular configuration is interference-free, and (2) planning a path between interference-free configurations.How and when these tasks are performed is determined by the particular path-planning strategy employed.These strategies can be categorized as either global or local, based on how they explore the configuration space.
Global path planning strategies first construct and partition the entire configuration space into cells, individually assessing each cell for interference.Path planning is performed only once the configuration space has been fully mapped out.These strategies are effective when the topography of the configuration space can be captured with a reasonably coarse sampling.However, for tight-fitting insertion problems, where the shapes are complex and the clearance is small, a prohibitively fine sampling may be required in the vicinity of the inserted configuration.This is exacerbated by the high computational cost of the methods used to assess interference, which typically involve dividing the geometry into fine polyhedra and repeatedly testing for overlap between the polyhedra of different bodies.Moreover, for tight-fitting bodies, where small incremental movements are required for insertion, simplification of the geometry or reduction of the configuration space is not possible.Therefore, techniques such as hierarchical decomposition (Faverjon 1984;Zhu and Latombe 1989), geometry simplification (Hwang et al. 2003), or selective exploration of configuration space (Kavraki and Latombe 1994;Karaman and Frazzoli 2011) are not applicable, meaning that complete mapping of the configuration space remains infeasible.
Local path planning strategies reduce the number of interference assessments by starting from an initial configuration and using local information to move toward the intended final configuration, mapping out the configuration space only as needed.Conventional local path-planning algorithms, such as artificial potential fields, Dijkstra, A*, D*, and rapidly exploring random trees (Siciliano et al. 2008), are often able to plan and adapt a path based on new environmental information but tend to be computationally expensive or become unstable in tight-fitting environments (Ayawli et al. 2018).Intelligent reactive approaches search the configuration space more efficiently by leveraging natureinspired principles (Mac et al. 2016).Among these, genetic algorithms, artificial neural networks, and fuzzy logic have proven effective in navigating some specific types of environments (Siddique and Adeli 2015).
Despite these advances, few path-planning strategies are efficient in navigating tight-fitting environments with complex geometry.One exception is a method proposed by Joskowicz and Taylor (1996), which specifically targets the peg-in-hole problem by exploiting the tight fit and known insertion direction.It works by partitioning the cavity into neighborhoods.These neighborhoods, defined by convex polyhedra, are then used to generate a set of linear constraints that ensure non-interference, thereby eliminating the need for overlap tests.The algorithm seeks out a sequence of small movements by formulating, for each successive movement, a linear optimization problem in which a move-related task function is minimized while subject to a set of neighborhood constraints.This method significantly reduces the cost of each movement.Yet, like other pathplanning approaches, it involves searching for a complete path, which may be made up of many individual configurations, each adding to the computational cost of the insertability assessment.
Related to insertability is the problem of ensuring accessibility for machining.Langelaar (2019) suggested a constraint developed for multi-axis machining.This method assesses whether a particular design can be machined from a block of material by effectively making sure that any material to be removed is accessible from one of several prescribed insertion directions.If a single insertion direction is considered, any cavity which can be created by the tool can be interpreted as a hole for which there exists a complementary peg.This method is highly efficient since the insertion path is known.Moreover, it provides gradient information that can be used in design optimization.Nonetheless, it is effective only when an exact rectilinear insertion direction is provided and rotation along the insertion path can be omitted.
In this paper, we present a novel approach that exploits the tight-fitting condition to assess insertability without computing the entire insertion path.We represent insertability as a continuously differentiable function that can be incorporated into gradient-based optimization to ensure that, to the extent that is possible, the optimized design satisfies the insertability requirement, together with the other design specifications.The insertability analysis can, then, be performed throughout the design process, rather than as a validation tool.

Page 3 of 19 220
We formulate the insertability problem as a path-planning problem in which the object is initially fully inserted.The objective is to determine if there exists a path from this inserted configuration to some uninserted configuration.We eliminate the need to construct the full path by exploiting the fact that the body and cavity shapes are nearly complementary, which is typically the case in the design and optimization of implants.We prove that if there exists any infinitesimal interference-free movement which moves the body slightly out of the cavity, then the object can be fully extracted, without interference.We demonstrate how our formulation can be employed as a design tool for both generic shape optimization and for advanced structural optimization.
The remainder of this paper is organized as follows.Section 2 describes the novel insertability metric and the related interference-free criterion.Section 3 presents a shape optimization strategy that incorporates the insertability metric as a design constraint.Section 4 explores an application of the insertability constraint to the optimization of a patientspecific acetabular cup for use in hip replacement surgery.Finally, Sects.5 and 6 conclude with a discussion and summary of the numerical results.

Insertability
Consider the perfectly complementary body-cavity system shown in Fig. 1.If the body is slightly rotated, as shown, clearance will be created in some areas (green) and interference will be introduced in others (red).In areas with clearance, points along the body surface are displaced in the direction of the local cavity surface normal, while in areas with interference, they are displaced in the negative normal direction.If all points along the body surface move in the direction of the local cavity wall normal, the configuration is interference-free.The question thus becomes whether there exists such a rigid body motion.This observation eliminates the need for the computationally expensive geometry subdivision and overlap testing commonly used to assess interference.The following describes how interference is assessed mathematically.
We define the movable object B as a rigid 3D body with 6 degrees of freedom (DOF), and the open cavity C as a fixed rigid 3D obstacle, as illustrated in Fig. 2. We further define the cavity interior I as the empty space within C, and the cavity exterior E as ℝ∖I .Finally, the body-cavity interface Γ includes the complementary surfaces of B and C, which are assumed perfectly coincident.We define a fixed global coordinate system and a movable frame associated with B.
We describe a configuration as the set (p, ) , made up of three displacements, p , and three rotations, , with respect to the fixed coordinate frame.
In a configuration (p, ) , the position, v , of a body point, b is expressed as where Rot( ) is the rotation matrix specifying the orienta- tion of the body.C o n s i d e r n o w t h e r i g i d b o d y m o t i o n T ∶ (p 0 , 0 ) → (p * , * ) , which maps from an initial configu- ration (p 0 , 0 ) to a new configuration (p * , * ) .The resulting displacement for a point b is: If point b is on the body-cavity interface, Γ , in (p 0 , 0 ) , r b may be expressed in terms of local interface normal and interface tangent components, ̂ b and ̂ b , as illustrated in Fig. 3: where nb is defined as pointing toward the cavity interior.

Interference-free criterion
Interference is defined as the condition in which two bodies overlap at some point in space.For a continuous motion, starting from a tight-fitting configuration, interference can only initiate where both bodies were previously coincident.Therefore, analysis of the instantaneous velocities of the bodies along their interface can predict subsequent interference.
Let us now consider a motion T ∶ (p 0 , 0 ) → (p * , * ) which results in an infinitesimally small movement of the (2) body.If the resulting displacement r b of a point b on the interface Γ is positive with respect to the local normal, i.e., then it can be said that T is collision-free with respect to b.For the general case in which the tangential displacement is non-zero, the displacements should be limited based on the shared local interface curvature, b .If the local curvature is less than or equal to 0, the minimum normal displacement is unrelated to the tangential displacement.However, when the local curvature is positive, a purely tangential displacement causes interference.In this case, a minimum normal displacement r n min is defined based on the local interface curvature, , and the tangential displacement r t b , as shown on the bottom right of Fig. 3. Graphically, r n min is the distance in the normal direction nb from v 0 b to the intersection of the local curvature circle and a line parallel to r n b through the point v 0 b + r t b tb .Mathematically, r n min is given by: and Eq. 4 is replaced by the general local Interference-free criterion (4) Since we are concerned with infinitesimal displacements, the local curvature is not expected to change significantly from v 0 to v * .The constraint is, therefore, sufficient to prevent interference.
If Eq. 6 is respected for all points along Γ , then T is inter- ference-free.This can be expressed as: To facilitate differentiability, the Kreisselmeier-Steinhauser (KS) function commonly employed in numerical optimization (Kreisselmeier and Steinhauser 1980;Akgun et al. 1999;Poon and Martins 2007) is used to approximate the maximum function.The final Interference-free criterion is then expressed as: where is a tuning parameter, typically in the range 100-1000.

Path planning
Traditionally, insertability analysis requires the iterative computation of a complete path between the inserted and extracted configurations.The high computational cost associated with this process is exacerbated when full 6 DOF motion and complex geometry are considered.Here, we exploit the tight-fitting condition to circumvent the process entirely.
Consider the perfectly complementary body-cavity system shown in Fig. 4 and an interference-free transformation T ∶ (p 0 , 0 ) → (p * , * ) .T is defined as non-reentrant if it results in an exclusively outward movement of the body with respect to the cavity.In other words, any point b on the cavity opening O in the initial configuration must either remain on O or move to the cavity exterior E, as a result of the movement T , i.e., If there exists a non-zero T which simultaneously satisfies the non-reentrance and the interference-free criteria, it can be proven that successive application of the same transformation T on B will always be interference-free and will eventually lead to full extraction from the cavity.In other words, T describes a full interference-free extraction/inser- tion path.

Property 1 Successive application of T on B is always interference-free
Proof Consider a body point b j which initially lies inside the cavity interior I.The transformation ), since T satisfies the non-reentrance criterion.Otherwise, if v * j ∈ I (case 2), then v * j is coincident with some position v k corresponding to a body point b k in the initial configuration (p 0 , 0 ) .Since T (p 0 , 0 ) → (p * , * ), b k was interference- free, T (p * , * ) → (p * * , * * ), b j is guaranteed to be interfer- ence-free.Thus, repeated application of T on B will always be interference-free.

Property 2 Successive application of T on B eventually leads to complete extraction
Proof Excluding the trivial zero-displacement solution, any transformation T which satisfies the non-reentrance constraint results in a portion of the body exiting the cavity.Suppose T This means that each application of T on B causes some portion of B to exit the cavity.Since B is always purely exiting, repeated application of T will eventually lead to the full extraction of B.
The process of assessing whether a body-cavity system is insertable along a path described by the movement T is described in Algorithm 1.The algorithm reports the degree of insertability, measured by the amount of interference (Eq.8.This scalar-valued metric allows an optimization routine to find the motion which minimizes interference.

Algorithm 1: Insertability
Input: A set of n interface vertex coordinates Γ x , unit normals Γ n, and curvature values Γ κ , an immersed condition function, ψ, and a transformation T Output: A scalar value I representing the degree of insertability of design x with respect to the transformation T

Design optimization
The proposed insertability metric, defined as the scalar output of Eq. 8, can be used as an efficient design validation tool.However, it is most effectively used as a constraint in structural optimization.To demonstrate this, we formulate a shape/topology optimization problem, following a common level set approach (Andreasen et al. 2020).This formulation is ideally suited to the interference constraint since it provides a crisp interface and allows for an efficient representation of arbitrarily complex topology (Van Dijk et al. 2013).
The following sections describe the level set method and provide detail on how the insertability constraint is implemented.

Design parametrization
The continuous body-cavity system is represented as a discrete nodal level set field on a regular grid, as shown in This representation provides a crisp interface and automatically generates a set of sample interface points for assessing insertability, along with their respective surface normals and curvature.

Interface properties
The zero contour provides a polygonal representation of the smooth geometry.In spite of this, the local normal and curvature information at the vertices can be readily obtained from the gradient of the level set field: (10) where l−1,m,n and l+1,m,n are the neighboring level set values in the x-direction.The other terms are calculated accordingly for each physical dimension.The nodal unit normal vectors are obtained by normalizing the local gradient, i.e., and the curvature is obtained as the divergence of the normalized gradient field: The zero contour unit normals n and curvature are interpolated from n and via Eq.10.

Partitioning
Finally, the body level set zero contour vertices are partitioned into a set of body-cavity interface vertices and a set of external vertices.This is achieved by converting the immersed region (cavity body and interior) into a level set (x) , and interpolating its value at each body contour vertex.Positive values correspond to interface vertices and negative values correspond to external vertices.
The discrete field of design variables is not used directly as a level set.Instead, we follow the filter-projection-scaling method presented in (Andreasen et al. 2020) to ensure stability and mesh independency.The three stages are described briefly here, for completeness.

Filtering
Firstly, the raw field of design variables, s , is filtered to s .The standard Helmholtz filter, introduced by Lazarov and Sigmund (2011), produces a smoothed field s as the solution of the equation with homogeneous Neumann boundary condition s n = 0 imposed on the boundary of the design domain.The radius parameter r defines the smoothing strength.

Projection
The projection stage applies the smoothed Heaviside function frequently used in robust topology optimization (Wang et al. 2011), to sharpen the filtered design.This stage helps stabilize the optimization procedure and can be used to generate additional eroded and dilated variants, in the case of robust optimization.
where is a threshold parameter, and controls the steepness of the projection.Throughout this work, = 0.5 and = 12.

Scaling
The final scaling step ensures that ∈ [ min ; max ] , through linear rescaling, i.e., The interval is chosen as ∈ [−h∕2;h∕2] , where h is the side length of a grid element, following Sharma and Maute (2018) and Andreasen et al. (2020).

Shape modification for insertability
As a first example, we consider a number of 2D and 3D bodycavity systems and use the proposed insertability metric to assess whether they are insertable and, if not, determine how they can be modified to ensure insertability.To that end, we formulate an optimization problem wherein the objective is to minimize a cost function, defined as total shape change (Eq.18, and the insertability requirement (Eq.8) is included as a constraint.This formulation ensures that the design will be rendered insertable while changing the shape as little as possible.The optimization routine then aims to simultaneously find the optimal shape, together with the optimal insertion motion, as defined by s and T.
The shape match function is defined as: where d(b i , S 0 ) is the closest distance from a point b i on the evolving body surface to a point cloud describing the original geometry S 0 , which can be interpolated from a pre-com- puted signed distance field.The resulting constrained optimization problem is solved by mathematical programming using the method of moving asymptotes (MMA) (Svanberg 1987).The procedure for optimizing a design for insertability is shown in Algorithm 2. It is also necessary to exclude the trivial solution with transformation T resulting in a triv- ial zero displacement, as it does not guarantee insertability.This is achieved by imposing a constraint on the minimum net displacement of the body subject to T.
Algorithm 2: Make insertable Input: A discrete design field s 0 , an immersed condition function ψ Output: A modified design field s In the first case, the initial design is axisymmetric and the interference is localized around the undercuts, as shown.The modified design eliminates the undercuts to produce a straight shaft.In the second case, doublecurvature is present, resulting in three regions with local interference.The algorithm simultaneously eliminates the double-curvature and modifies the extraction direction to eliminate the interference.The third case is similar to the first and is included to show that the algorithm can handle complex topology such as internal cavities.In this instance, a flooding algorithm (Heckbert 1990) is used to distinguish the internal and external contours.Finally, case four highlights the algorithm's ability to handle large rotations.In this instance, the contour is modified to smooth out the protrusion, while largely maintaining its curvature.

Space-filling
Another noteworthy feature of these results is that the algorithm tends to move interface inwards in some areas and outwards in others.In cases where both the body and cavity can be modified, this is adequate.However, in many cases, the cavity geometry is prescribed, and the objective may be to design a body that maximally fills the cavity.In this case, the space-filling version of the shape match function may be used: where and are tuning parameters, with  >>  .This ensures that outwards changes are more expensive than inward changes.The same original designs presented in Fig. 6 are modified using the shape-filling objective.Results are shown in the top row of Fig. 7.In each case, the cavity is unchanged and the design is reduced in order to comply with the insertability constraint.

Space-making
We also consider a scenario in which only the cavity can be modified, or in which the body must at least occupy the same space as the original.The direction-distinguishing version of the shape-matching function Eq. 19 can be used with  <<  .Results obtained using this space-making version are presented in the bottom row of Fig. 7. Here, the results refer to the scenario where both body and cavity can be modified.In the case where only the cavity can be modified, the modified interface would be as shown, but the body would remain unchanged.

3D design
The same method is equally applicable to the modification of 3D designs.Various 3D original and modified designs are presented in Fig. 8.In 3D, more complex movements are possible.For example, the third example involves a corkscrew motion about the vertical axis.The algorithm successfully identifies the appropriate transformation vector automatically.As for the 2D case, the space-filling and space-making versions may be used instead of the original shape-matching function.Space-making and space-filling results are presented in the lower half of Fig. 8. Numerical results for all studies are included in Table 1.

Ease of insertion
The insertability constraint is defined as an inequality with the right side set to zero.This, in effect, allows for pure sliding motion along the body-cavity interface.In certain applications, such as in automated assembly, this can result in designs that are difficult to insert.If, for example, the body and cavity have rough surfaces, the resulting friction may require high forces or prevent insertion entirely.This can be addressed by setting the right side of the inequality constraint to ≤ 0 .Figure 9 shows how the design of a straight bar is affected for various values of .

Filter radius
A minor issue highlighted by these results is the unnecessary dilation near the cavity opening.This results from the filtering stage described above.While it may be possible to eliminate this issue by increasing the tuning parameters in the projection stage, it is generally sufficient to replace the exterior portion of the modified design with that of the original, as shown in Fig. 9d.We note that in some cases, the modified design does not maximally fill the cavity.In cases 2 and 3, for example, the implant should touch the bottom of the cavity, but does not.This is likely related to the filtering and may be addressed by either adjusting the tuning parameters or refining the levelset grid.

Convergence
The convergence history for a selection of the case studies, chosen at random, is shown in Fig. 10.In all cases, the interference constraint is initially normalized.The optimization history reveals a generally smooth evolution in which the design is first changed significantly at the expense of the objective function, in order to bring the interference constraint below zero.Subsequent iterations minimize the shape match objective.The optimization tends to converge within less than 100 iterations.Note that in the 2D case (Fig. 7a), the observed instability between iterations 75 and 110 results from a coarse mesh discretization and a large step size.

Application to the design of acetabular cups
Orthopedic implants are medical devices used to replace damaged or missing bone tissue in an effort to improve the patient's mobility, reduce pain, and generally improve quality of life.Traditional implant systems include a set of generic implant shapes and sizes to help find the best fit.However, patient-specific implants (PSI), designed according to a patient's unique anatomy, may offer better outcomes than their traditional counterparts.Some of the expected benefits of patient-specific implants include: better positional accuracy of the joint, reduced risk of complications, and reduced recovery times (Small et al. 2014;Spencer-Gardner et al. 2016;Maniar and Singhi 2014).
A particularly challenging aspect in the design of patientspecific implants is the need to optimize functional compatibility between the implant and the surrounding bone tissue, while also ensuring that the implant can be inserted without causing damage to itself or to the peri-prosthetic tissue.The method presented in this work is ideally suited to address these potentially conflicting design objectives.In this section, we demonstrate how our insertability constraint can be applied to the design of an acetabular cup, considering both geometric and structural objectives, via finite element analysis (FEM).
The acetabular cup is the female component of a hip replacement implant (Fig. 11).It is a typically hemispherical component that is placed in the acetabulum, which is the socket of the pelvis that forms the hip joint.Traditionally, orthopedic surgeons determine the ideal diameter of the cup from a range of options, based on patient-specific measurements.However, since the surface of the acetabulum is not perfectly hemispherical, there is always some morphological mismatch between the cup and the acetabulum.This mismatch may reduce the implant's stability and cause damage to the surrounding cortical bone.
As a first approach, we can apply the shape-matching formulation subject to our insertability constraint, as described above.However, we must first generate a morphologically matched implant design to use as a starting point.This can be achieved by performing a CT scan of both hips and extracting the difference between the healthy and defective anatomy.If both hips are anatomically defective, we may use a computational model to obtain an approximate reconstruction of the original healthy bone.In this work, we use a pre-processed set of defective and reconstructed pelvis geometry obtained from Meynen et al. (2020).The defective pelvis, which is classified as type IIA according to the Paprosky classification system, is reconstructed using a statistical shape model based on data from 90 patients.By comparing the defective geometry, shown in Fig. 12a with the reconstructed healthy geometry, shown in Fig. 12b, we can identify the missing acetabular tissue and the ideal bone-implant contact surface.The missing tissue geometry is then used to create an implant with a perfectly matching bone contact surface, as well as the required hemispherical bearing surface for the synthetic joint, as shown in Fig. 13.
Results of the design process are presented in Fig. 14.Each column shows a different view of the implant-bone contact surface.The top row shows the values of the local insertability metric.The areas which prevent insertability are highlighted with circles.The bottom row shows the optimized implant with the local shape change shown in red and green.In this case, since cortical bone tissue can be resected but not added, green areas represent gaps between the implant and the bone.We can observe that the shape changes are concentrated on either side of the implant and that, for each, inwards and outwards changes are paired.This coupling makes the local surface parallel to the insertion direction, as would be expected for the nearly linear insertion path shown.
If indicated, space-filling and space-making may be used instead.The space-filling design will preserve as much of the cortical bone tissue as possible but will result in some morphological mismatch.On the other hand, the space-making design will ensure a perfect contact surface but requires some bone resection.
Beyond the purely geometric shape-matching functions, more sophisticated objectives based on structural analysis of the bone-implant system may be considered.For example, we may aim to minimize the risk of implant slipping and loosening by assessing the global Hoffman failure criterion Hoffman (1967).This function assesses the local risk of interface fracture based on the cortical bone properties and the local stress conditions throughout the interface.With the Hoffman index set as the optimization objective and the insertability criterion set as the only constraint, and reduction in the Hoffman index of 94% is achieved (Fig. 15).The Hoffman criterion and the associated finite element methods

Discussion
The preceding examples demonstrate that the proposed insertability constraint can effectively and efficiently restrict the design of complex 2D and 3D structures subject to various geometry and structural requirements.
The method is implemented in MATLAB R2023a.All 2D and 3D results presented were computed on a laptop computer with M1 Max processor and 64 GB LPDDR5 RAM.The 2D shape optimization was performed with approximately 20,000 design variables and required 30-50 iterations, each performed in approximately 0.6 s.The 3D results were performed with approximately 500,000 variables and required 40-100 iterations, at 8-10 s per iteration.Compared with our MATLAB implementation of the method presented by Joskowicz and Taylor (1996), the 2D insertability assessment was performed approximately 70× faster, primarily due to the computation of a complete insertion path.
The level set-based optimization framework proved highly effective as a host for the insertability constraint.Performance improved monotonically throughout the optimization procedure and converged in less than 100 iterations for all the experiments.The algorithm proved capable of handling complex topology and cases with large rotations in both 2D and 3D.No parameter tuning was required between numerical studies.

Limitations
The main limitation of the insertability metric is that it reduces the cavity to an interface curve (2D) or surface (3D), thereby ignoring possible interference between the fixed cavity body and the external/extracted portion of the moving body.For this reason, under certain circumstances, the insertability criterion may potentially produce false positive results.Handling of this special case may require additional constraints and is left as a challenge for future research.

Notes on implementation
For problems in which significant rotation is required for insertion/extraction, it is critical to choose an appropriate reference frame origin.Otherwise, the sensitivities of the constraint with respect to the rotations in may become very large.This is generally avoided by setting the origin to the centroid of the initial design.
The insertability constraint tends to be sensitive to the filter radius.In particular, significant instability has been observed for r < 4 .While a relatively large filter radius does tend to reduce the design complexity for a given grid size, sharper projection and finer grid resolution help mitigate this issue.
In a real-world application, designs should take into account manufacturing-related uncertainties.Methods such as the one presented by da Silva et al. (2019) help account for local defects on internal stresses.In the case of the proposed insertability constraint, small surface imperfections may impede insertion.We demonstrated how this can potentially be mitigated by tailoring the ease of insertion parameter .However, a perfect fit in the inserted configuration remains tied to the accuracy of the manufactured components.For this reason, manufacturing uncertainty should be incorporated into the optimization framework, alongside the insertability constraint.

Conclusion
We have demonstrated that the proposed insertability metric is able to efficiently assess the insertability of complex body-cavity systems, taking into account full 6 DOF motion.Beyond its use as a validation tool, the insertability metric can be applied as a constraint in the design and optimization of functional components where geometric or structural objectives are considered.The level set-based formulation allows for straightforward integration with many topology optimization approaches.Future work will focus on applying the insertability constraint to multi-scale optimization in which both the overall shape and the internal microarchitecture of the architected material are optimized simultaneously.

Appendix A: Finite element methods
The finite element analysis used to assess the risk of interface fracture is based on an enriched XFEM model with second-order displacement discontinuities along the boneimplant interface.The fixed isometric grid is extended from the design level set grid to encompass the entire pelvic bone and is comprised of tri-linear hexahedral elements with 1.45 mm side length.A detailed description of the XFEM strategy can be found in Moës et al. (1999).
In order to assess the stress conditions along the boneimplant interface, appropriate loads and boundary conditions are essential.Ghosh et al. (2015) proposed modeling both the sacroiliac joint and pubic symphysis as fully constrained.Instead, we apply an in-plane constraint on the pubic symphysis, which, we argue, more accurately reflects the bilateral symmetry of the full pelvis.The joint reaction force is distributed on the implant's inner surface, following the elastic body bearing pressure model described in Popov et al. (2019).All loads and boundary conditions are shown in Fig. 16.The total magnitude and orientation correspond to the walking case described in Heller et al. (2005).
Approximate bone densities and mechanical properties are mapped from the CT scan.Densities are obtained by linearly rescaling the Hounsfield values to a typical range of [0.04; 1.9] g cm −3 in accordance with Dalstra et al. (1993).The elastic moduli are then interpolated from the densities according to Morgan et al. (2018).The densitydependent material strengths used to assess the Hoffman failure risk index are computed as per Stone et al. (1983).

Appendix B: Sensitivity analysis
In order to use the insertability function in an optimization setting, we require the sensitivities with respect to the field of design variables.They are computed analytically with respect to the processed design variables as: where the partial sensitivities of the local interference function at each vertex i are:  Note that these sensitivities are computed with respect to the processed design variables .To obtain the sensitivities with respect to the actual design variables, the chain rule is applied.

Fig. 1
Fig. 1 A body-cavity system shown before and after a small body movement.Clearance and interference are shown in green and red, respectively.Local cavity surface normals are shown in yellow.Sample body surface vertices in areas with interference are displaced in the negative normal direction, while body surface vertices in areas with clearance are displaced in the positive normal direction.(Color figure online) Fig. 2 Illustration of an original (left) and partitioned (right) bodycavity system

Fig. 3
Fig.3A body-cavity system in its inserted configuration, C 0 , (left) and a perturbed configuration, C * (center).The net displacement of a point b from v 0 to v * of an arbitrary interface point (red) is decomposed into

Fig. 4
Fig. 4 An illustration of both possible scenarios resulting from the application of transformation T to a body point b j .In case 1, b j moves from v 0 j to v * j , which lies outside of the cavity interior.In case 2, b j moves from v 0 j to v * j , which overlaps with the position v 0 k of some point b k in the initial configuration C 0 Fig. 5.The surface geometry is defined as the zero contour of the level set , with  i > 0 designating the interior of the body, and  i < 0 designating the exterior.The positions of the zero contour vertices x are linearly interpolated from the nodal level set values according to:

Fig. 5
Fig. 5 An illustration of a continuous body-cavity system (left), approximated as a discrete level set (right).The level set zero contour is partitioned into interface vertices (red) and external vertices (black).By convention, > 0 on the interior and < 0 on the exterior.(Color figure online) Γκ,ψ,T ) 13 c d ← −DISP LACEMENT (Γx,T ) 14 o s ←SHAP ECHANGE(Γx,Γ 0 x ) 15 (s, T ) ←MMA(min (s,T ) (os), s.t.{ci, c d }≤0) 16 δ ← |(s, T ) − (s, T ) prev | 17 while c i > 0 | c d > 0 | δ > δ min 18 return s Four uninsertable designs are shown in the top row of Fig. 6.The corresponding optimized designs are shown in the bottom row.In each case, the extraction direction is shown by an arrow, and the areas which locally violate the insertability constraint are highlighted in red.Once the iterative optimization procedure has converged, the shape has been modified to eliminate all local interference.Magenta lines highlight the local shape change from the original design.

Fig. 6
Fig. 6 Various 2D body-cavity systems modified for insertability using Eq.18 as the design objective.The original designs (top row) fail to meet the local interference-free criterion at the highlighted points (red).The modified designs (bottom row) are presented along with the local shape change from the original design, shown as magenta vectors.The crosshatching designates the cavity wall.All designs use 2 × 10 4 node level set meshes with a filter radius of r = 6 .(Color figure online)

Fig. 7 6 220
Fig. 7 Results obtained by applying the space-filling (top row) and space-making (bottom row) objectives to the original designs from Fig. 6

Fig. 8
Fig. 8 Various 3D body-cavity systems modified for insertability using Eq.18 as design objective.The original designs (top row) fail to meet the local interference-free criterion at the highlighted areas (left).The modified designs (bottom rows) are presented along with the local shape change from the original design (right).Red and green

Fig. 9
Fig. 9 a-c Three versions of a straight rod, each modified for insertability with different values.From left to right, = 0 , = −1 × 10 −2 , = −1 × 10 −1 .d Shows how the exterior portion of the design in c can be overwritten with the original geometry from a

Fig. 11
Fig. 11 An illustration of the total hip replacement (THR) procedure with a generic hemispherical acetabular cup.Reaming ensures proper contact between the bone and cup along the exterior, but pockets are

Fig. 12 a
Fig. 12 a A Paprosky classification type IIA defective pelvic bone extracted from CT data, b the reconstructed healthy bone obtained using the statistical shape modeling by Meynen et al. (2020, c the defective pelvic bone with the recovered acetabular tissue highlighted

Fig. 13
Fig. 13 The initial implant level set based on the geometry of the missing acetabular tissue.The bone-implant interface is highlighted in blue while the modified spherical bearing surface, required for finite element analysis, is highlighted in red.(Color figure online)

Fig. 15
Fig.15The original (top row) and modified (middle row) implant designs colored according to the local risk of interface fracture fromGarner et al. (2022).In the bottom row, the local shape change is sensitivities of the interpolated vertex position b i and gradient ∇ i are equal to zero everywhere except at their respective edge parent nodes p1 and p2, for which: and where x represents either b i or ∇ i .

Fig. 16
Fig. 16An illustration of the loads and boundary conditions applied on the pelvis-implant system.The pelvis is fully fixated along the sacroiliac joint, and an in-plane constraint on the pubic symphysis prevents motion in the lateral direction.The implant-acetabulum interface is modeled as fully bonded.A joint reaction force corresponding to the walking case described inHeller et al. (2005) is distributed on the spherical contact surface, following the elastic body bearing pressure model described inPopov et al. (2019)

Table 1
Problem parameters and performance values for each numerical experimentThe insertability metric values are normalized with respect to the mean body displacement, and the Hoffman risk index is normalized to that of the original design