# Co-generation of Collision-Free Shapes for Arbitrary One-Parametric Motion

Clinton B. Morris and Morad Behandish  Palo Alto Research Center (PARC), 3333
Coyote Hill Road, Palo Alto, California 94304

###### Abstract

Mechanical assemblies can exhibit complex relative motions, during which
collisions between moving parts and their surroundings must be avoided. To
define feasible design spaces for each part’s shape, “maximal” collision-free
pointsets can be computed using configuration space modeling techniques such
as Minkowski operations and sweep/unsweep. For example, for a pair of parts
undergoing a given relative motion, to make the problem well-posed, the
geometry of one part (chosen arbitrarily) must be fixed to compute the maximal
shape of the other part by an unsweep operation. Making such arbitrary choices
in a multi-component assembly can place unnecessary restrictions on the design
space. A broader family of collision-free pairs of parts can be explored, if
fixing the geometry of a component is not required. In this paper, we
formalize this family of collision-free shapes and introduce a generic method
for generating a broad subset of them. Our procedure, which is an extension of
the unsweep, allows for co-generation of a pair of geometries which are
modified incrementally and simultaneously to avoid collision. We demonstrate
the effectiveness and scalability of our procedure in both 2D and 3D by
generating a variety of collision-free shapes. Notably, we show that our
approach can automatically generate freeform cam and follower profiles, gear
teeth, and screw threads, starting from colliding blocks of materials, solely
from a specification of relative motion and without the use of any feature-
informed heuristics. Moreover, our approach provides continuous measures of
collision that can be incorporated into standard gradient-descent design
optimization, allowing for simultaneous collision-free and physics-informed
co-design of mechanical parts for assembly.

###### keywords:

Automated Design , Collision Avoidance , Configuration Space , Spatial
Reasoning , Persistent Contact

††journal: Computer-Aided Design

##  1 Introduction

Design for manufacturing and assembly (DfM/DfA) [1], popularized by Boothroyd
and Dewhurst in the early 1980s, is a set of principles for design of
manufacturable parts to enable faster, cheaper, and more reliable assemblies.
These principles include minimizing the number of parts and assembly
operations, making these operations simple and fail-proof, standardizing the
interfaces, etc. Such design rules and evaluation criteria can be effective as
guidelines for manual design or fairly limited automation with parameterized
geometric features that can be manufactured with traditional processes. With
the advent of modern (e.g., additive and hybrid [2]) manufacturing, there are
tremendous opportunities for DfM/DfA demanding a deeper look into geometric
and physical modeling of assemblies.

Computational design tools such as gradient-based shape and topology
optimization (SO/TO) [3, 4] have proven effective to satisfy physics-based
performance criteria at the individual part level, owing to their formulation
as differentiable objective functions and constraints. However, incorporating
assembly and manufacturing constraints, often formulated using geometric and
kinematic models, is challenging [5]. A quantitative analysis of concepts such
as collision, containment, contact, and complementarity of parts of
arbitrarily shapes undergoing arbitrary motions is required, before they can
be cast into differentiable measures. This paper focuses on collision and
containment, quantified by overlap measures between shape indicator functions
[6], and sets the groundwork of an extensible framework to more complex
concepts such as contact and complementarity, quantified by more complex shape
functions [7, 8].

Consider, for example, the case of a single degree-of-freedom (DOF) mechanism
such as a four-bar or slider-crank linkage, a cam/follower or pinion/gear
pair, a latching mechanism, etc. The design of such assemblies often starts
from a system-level design (e.g., determining types and relative positions of
joints on each link) to satisfy kinematic and dynamic requirements,
determining the relative motions between parts [9]. While these computations
are automated and available in commercial packages, there is a gap in
cascading them to part-scale 3D design with any level of generality. For
example, there is no systematic way to impose collision avoidance between the
crank and coupler in a four-bar linkage or persistent contact between the cam
and follower with prescribed motion, while optimizing their shapes for
performance (e.g., stiffness and strength) subject to manufacturing
constraints [10, 11, 12]. Fig. 1 [13] illustrates a more complex assembly with
various collision, containment, contact, and complementarity requirements that
are difficult to formalize for design automation. Ultimately, to automate the
design or redesign of such a complex assembly, all of these concepts must be
integrated into automated design workflows. While designing such a complex
assembly is beyond current capabilities, automated collision avoidance is a
necessary advancement to achieve this goal.

Figure 1: Cross-section of engine demonstrating collision, contact,
containment, and complementarity constraints to properly function [13].

In this paper, we devote specific attention to collision avoidance or
(equivalently) containment constraints, as they capture a broad class of
design criteria involving interactions of shapes and motions [14, 15]. To
incorporate such constraints into gradient-based optimization, the extent of
their violations must be locally evaluable and differentiable with respect to
design variables (e.g., small geometric or topological changes to the design),
to penalize the objective functions. Additionally, collision is a pairwise
relation, presenting additional challenges when multiple moving parts are to
be designed simultaneously. Changes to one part will potentially introduce or
eliminate collisions with the others, necessitating a co-design workflow in
which the parts are shaped simultaneously, evolution of one part directly
informing those of the others in the assembly.

To the best of our knowledge, such a co-design workflow that enables designing
complex shapes under arbitrary motions does not exist. This paper presents a
framework to develop such workflows to design families of parts in which
collision avoidance constraints can be simultaneously satisfied and seamlessly
integrated with other (e.g., performance and manufacturing) constraints. More
specifically, we present a procedure to generate members of a family of
“maximal pairs” of collision-free parts. We do not divert our attention to
formulating other constraints in this paper, although we show a general recipe
to couple them with collision constraints.

###  1.1 Related Work

Existing approaches to generating shapes that satisfy motion-based collision
and contact constraints place excessive restrictions on the design space,
missing opportunities for better-performing and more cost-effective
assemblies. Ilieş and Shapiro [14, 15] developed a framework to produce a
“maximal” shape for a part, moving against another part of prescribed shape,
to satisfy collision avoidance or (equivalently) containment constraints. This
led to the definition of a fundamental new solid modeling operation called
unsweep [16, 17]. Nelaturi and Shapiro [18] extended the idea to a broader
class of configuration space operations, based on group morphology [19]. These
operations have proven effective in solving various manufacturing analysis and
process planning problems [20, 21, 22].

There are at least two challenges with using such operations, formulated in a
set-theoretic language, in a design framework. First, using them to compute a
maximal entity (shape or motion) requires full knowledge of all of the other
entities—for instance, sweep/unsweep map a given shape and a given motion to a
maximal shape, while homogeneous Minkowski products/quotients map a given pair
of shapes to a maximal motion, avoiding collisions [19]. Second, these set-
theoretic operations do not interoperate well with other (e.g., performance
and manufacturing) constraints, commonly expressed using analytic
inequalities. Unlike the latter, the former do not provide quantitative
measures of violation of collision constraints, which are critical for
navigating tradeoffs with other constraints. Measure-theoretic generalizations
of such operations [23] based on convolution algebras [19] can be used to
obtain locally evaluable and differentiable measures that can be integrated
with other constraints via Lagrange multipliers [5]. However, the attention
has so far been restricted to single part design.

In contrast, our procedure does not require either part’s shape to be fixed
upfront and can generate families of maximal collision-free designs (in a
partial ordering of pairs), using local and global measures that can be
integrated with other constraints and differentiated for gradient-descent
optimization. We present an incremental co-generation procedure that subsumes
one-way unsweeps of either part, against the initial design of the other
part’s complement, as extreme cases of the family of maximal pairs.

To avoid premature decisions and excessive geometric constraints on part
geometries, Stockli and Shea [24] developed an automated procedure for
simultaneously modifying the shapes of two parts in relative motion. Co-
generation of the shapes was achieved using a TO procedure in which they
formulated and minimized an aggregate collision measure. They showed that such
measures can be rapidly computed through matrix multiplications during the
optimization loop to globally quantify collisions. To update the shapes within
an iterative optimization loop to alleviate the collisions, basic rules were
applied to grow and shrink the shapes.

Our computational algorithm has similar elements to [24] in the way collision
measures are computed using a voxelization of the moving domains, although we
present a formal set-theoretic groundwork, built on top of [14, 15, 17, 18]
that is not tied to a specific representation scheme and can be discretized in
many different ways. Moreover, we use gradients and local measures of the
collision measures to augment the TO sensitivity fields, enabling further
scalability. Most importantly, our approach generates families of maximal
shapes as opposed to a single arbitrary pair obtained using rule-based
heuristics in [24].

###  1.2 Contributions & Outline

This article presents a general formulation and computational framework for
co-generating maximal collision-free shapes in arbitrary relative motion. We
show that:

  1. 1.

The unsweep operation and the underlying notion of maximality (for a partial
ordering of solids) can be expanded to pairs of solids to more broadly explore
the space of collision-free designs (Section 2).

  2. 2.

The collision of solids in relative motion can be measured locally and
globally to use as a differentiable violation measure to penalize gradient-
based optimization (Section 3).

  3. 3.

The violation measure can be used to formulate an iterative and incremental
co-generation procedure for a broad subset of maximal pairs of collision-free
solids (Section 4).

  4. 4.

The procedure can scale to efficiently co-generate a variety of nontrivial
collision-free solids that, if allowed by motions, exhibit persistent
contact111Persistent contact means nonempty boundary intersection maintained
throughout the motion. Here, contact is not explicitly enforced, although it
could be (e.g., by adding contact or complementarity measures to the objective
function). Its persistence cannot be guaranteed by shape design, as it also
depends on the given motion. When it is possible, maximizing volume subject to
collision-avoidance appears to improve its chances. Understanding persistent
contact properties of shapes and motions requires further research (out of
scope here). in 2D and 3D (Section 5).

  5. 5.

The violation measures enable using collision-avoidance constraints in
gradient-based optimization to extend beyond geometric reasoning (Section 6).

We also show that the collision measure can be efficiently computed by linear-
algebraic operations, involving a precomputed pairwise correlation matrix that
depends on the motion, pre- and post-multiplied by arrays of design variables,
all of which can be parallelized on the CPU/GPU (Section 3). However, this
separability of computations was first observed by Stockli and Shae [24],
hence does not constitute a novel contribution of this paper.

##  2 Preliminaries

Kinematic design of mechanical parts and assemblies under prescribed motions
can often be formulated as computing “maximal” pointsets that satisfy
collision avoidance or (equivalently) containment constraints. The maximality
is defined in the partial ordering of pointsets with respect to containment.
Such maximal shapes prune the design space to a feasible subspace, defined by
the powerset of the maximal pointset [5], for subsequent design space
exploration to satisfy various performance and manufacturing criteria (e.g.
via TO [10, 5, 11, 12, 25]).

Let $S\subset{\mathds{R}}^{\mathrm{d}}$ be a $\mathrm{d}-$dimensional solid or
‘r-set’, defined as a compact (bounded and closed) regular and semianalytic
pointset in the Euclidean $\mathrm{d}-$space ${\mathds{R}}^{\mathrm{d}}$ [26].
Let $M\subset{\mathrm{SE}(\mathrm{d})}$ be a motion acting on
${\mathds{R}}^{\mathrm{d}}$, i.e., a parameterized collection of rigid
configurations or ‘poses’ (combined rotations and translations) that a
$\mathrm{d}-$dimensional pointset can assume [27]. For computational purposes,
the solids are commonly represented by boundary or volumetric representations
(e.g., NURBS, CSG, surface/volume mesh, voxels, or sampled point clouds) [26],
while motions are typically represented by sampled or parameterized
homogeneous matrices, dual quaternions, etc. We restrict our attention to
$\mathrm{d}:=2,3$ and one-parameteric motions.

The first class of problems can be formulated and solved using the sweep and
unsweep operations [17]. The sweep of a given solid $S$ under a given motion
$M$ is another pointset $\mathsf{sweep}(M,S)\subset{\mathds{R}}^{\mathrm{d}}$,
which is a superset of $S$ collecting all points that are included in the
displaced shape $\tau S$ for at least one configurations $\tau\in M$:

  * •

Explicit definition (by an indexed union):

$\mathsf{sweep}(M,S)\triangleq\bigcup_{\tau\in M}\tau S.$

  * •

Implicit definition (by a membership test):

$\displaystyle\mathsf{sweep}(M,S)$
$\displaystyle\triangleq\big{\\{}{\mathbf{x}}\in{\mathds{R}}^{\mathrm{d}}~{}|~{}\exists\tau\in
M:{\mathbf{x}}\in\tau S\big{\\}}$
$\displaystyle=\big{\\{}{\mathbf{x}}\in{\mathds{R}}^{\mathrm{d}}~{}|~{}M^{-1}{\mathbf{x}}\cap
S\neq\emptyset\big{\\}},$

where $M{\mathbf{x}}\triangleq\\{\tau{\mathbf{x}}~{}|~{}\tau\in M\\}$ and
$\tau S\triangleq\\{\tau{\mathbf{x}}~{}|~{}{\mathbf{x}}\in S\\}$ where
$\tau{\mathbf{x}}\in{\mathds{R}}^{\mathrm{d}}$ denotes a displaced point.

Note that the membership of a query point
${\mathbf{x}}\in{\mathds{R}}^{\mathrm{d}}$ in the sweep is tested by applying
the inverse motion to the query point and checking whether the resulting
trajectory passes through the given shape.

On the other hand, the unsweep of a given solid $S$ under a given motion $M$
is another pointset $\mathsf{unsweep}(M,S)\subset{\mathds{R}}^{\mathrm{d}}$,
defined by the subset of $S$ including all points that that remain included in
the displaced shape $\tau S$ for all configurations $\tau\in M$:

  * •

Explicit definition (by an indexed intersection):

$\mathsf{unsweep}(M,S)\triangleq\\!\\!\bigcap_{\tau\in M^{-1}}\\!\\!\tau S.$

  * •

Implicit definition (by a membership test):

$\displaystyle\mathsf{unsweep}(M,S)$
$\displaystyle\triangleq\big{\\{}{\mathbf{x}}\in{\mathds{R}}^{\mathrm{d}}~{}|~{}\forall\tau\in
M^{-1}:{\mathbf{x}}\in\tau S\big{\\}}$
$\displaystyle=\big{\\{}{\mathbf{x}}\in{\mathds{R}}^{\mathrm{d}}~{}|~{}M{\mathbf{x}}\cap\overline{S}=\emptyset\big{\\}},$

where $\overline{S}\triangleq({\mathds{R}}^{\mathrm{d}}-S)$ is the set
complement of $S$.

Note that the membership of a query point
${\mathbf{x}}\in{\mathds{R}}^{\mathrm{d}}$ in the unsweep is tested by
applying the forward motion to the query point and checking whether the
resulting trajectory remains inside the given shape.

To see how sweep and unsweep can be used in practical design problems,
consider a pair of solids $S_{1},S_{2}\subset{\mathds{R}}^{\mathrm{d}}$ moving
under one-parametric motions $M_{1},M_{2}\subset{\mathrm{SE}(\mathrm{d})}$.
The solids can be parts in a single DOF mechanism, e.g., any pair of links in
a four-bar linkage, a cam/follower pair, a pinion/gear pair, a latch/pin pair,
etc. The relative motion (of $S_{2}$ as observed from a frame attached to
$S_{1}$) is $M\triangleq M_{1}^{-1}M_{2}$. If we fix the shape of solid
$S_{1}$, the maximal pointset that does not collide with $S_{1}$ (i.e., is
contained within $\overline{S}_{1}$) throughout the motion is:

$S_{2}^{\ast}\triangleq\overline{\mathsf{sweep}(M^{-1},S_{1})}=\mathsf{unsweep}(M,\overline{S}_{1}).$
(1)

The maximality implies that every other solid that satisfies the no-collision
constraint will be a subset of the maximal pointset (i.e., $(MS_{2}\cap
S_{1})=\emptyset$ iff $S_{2}\subseteq S_{2}^{\ast}$).

Let us assume the motion can be parameterized as $M=\\{\tau(t)~{}|~{}0\leq
t\leq 1\\}$. For the $1-$DOF mechanisms exemplified above, $t\in[0,1]$ can be
thought of as normalized time—or any parameter monotonically changing with
time—over the mechanism’s motion cycle. Consider the car hood latch example in
Fig. 2, revisited from [28].

Figure 2: Car hood latch example adopted from [28], in which a rotating latch
must remain contained in an enclosure while maintaining contact (without
collision) over a functional surface (in yellow) with a pin translating
vertically. The relative motion is shown on the left.

  * •

Our goal is to design the shape of a latch $S_{\text{L}}$, rotating around a
pivot by a given trajectory $M_{\text{L}}$ with respect to the car (common
frame of reference).

  * •

Constraint 1: Let $S_{\text{B}}$ represent a bounding region that the latch
must not exit to avoid collision with surrounding objects under the hood that
are stationary with respect to the car (hence $M_{\text{B}}$ is identity).

  * •

Constraint 2: Let $S_{\text{P}}$ represent a pin of given shape, attached to
the hood, thus moving along a given trajectory $M_{\text{P}}$ with respect to
the car. The latch must avoid collision while maintaining contact with the
pin, as it moves with a relative motion $M_{\text{P}}^{-1}M_{\text{L}}$.

The design space for the latch $S_{\text{L}}$ is thus restricted (by the
conjunction of these constraints) to subsets of

$S_{\text{L}}^{\ast}:=\mathsf{unsweep}(M_{\text{B}}^{-1}M_{\text{L}},S_{\text{B}})\cap\mathsf{unsweep}(M_{\text{P}}^{-1}M_{\text{L}},S_{\text{P}}).$
(2)

A similar approach can be used to design the cam profiles for given follower
profiles and cam/follower motions. Despite its generality (for arbitrarily
complex shapes and motions) and rigor, this approach has three major
limitations that we aim to address:

First, using unsweep requires pre-defining one of the two shapes in relative
motion to obtain a maximal shape for the other one. This requirement may be
fine for design problems with clear choices for the shapes of certain
components or spatial regions because of non-negotiable conditions (e.g., the
bounding box in Fig. 2) or simple profiles (e.g., the cylindrical pin in Fig.
2 and cylindrical/flat followers). However, in problems where both shapes in
relative motion are to be designed simultaneously (e.g., the cam/follower
system), the arbitrary choice of either shape limits the design freedom and
requires trial-and-error iterations over a theoretically uncountable—and
computationally exponential—number of possibilities.

Second, the maximal pointsets that satisfy these constraints come with no
guarantees for effective contact, shape complementary, force/torque
transmission, and other requirements, whose quantification for arbitrary
shapes is nontrivial. For example, consider two parts in a relative screw
motion, each to be selected as a subset of the bounding boxes. If one of the
parts is shaped as a cylindrical hole or a threaded nut, the maximal shape for
the other part will be computed by an unsweep as a matching cylindrical peg or
a threaded bolt, respectively. The contact and complementary is perfect in
this case, because the screw motions form a symmetry subgroup of
${\mathrm{SE}(3)}$ and the first part’s shape (the nut) was carefully selected
such that this particular screw motion of the same pitch is its symmetry;
hence unsweep generates a matching shape (the bolt) that, not only avoids
collision with the nut, but is also a perfect complement of it and remains as
such throughout the motion.222To be accurate, these statements apply strictly
to an infinite extension of a bolt/nut of the same pitch, meaning that
applying the motion to their shapes does not change the shapes (symmetry) and
they partition the $3-$space (complementarity) throughout the motion. However,
if we picked a slightly different shape for the first part, such as a half-
cylindrical hole or a nut of a different pitch, the unsweep would generate a
maximal shape that, despite avoiding collision, offers poor contact and
complementarity. This example illustrates why co-design is a more practical
approach where both shapes are chosen and refined simultaneously.

The third limitation is that set-theoretic operations such as unsweep, used to
express kinematic constraints, do not interoperate well with analytic
constraints such as global or local inequalities in terms of field variables,
used to express physics-based constraints (e.g., stiffness or strength). As a
result, concurrently satisfying kinematic and physical constraints is
challenging [22]. One possible approach is to translate set-theoretic
operations with measure-theoretic parallels, e.g., by replacing Boolean and
morphological operations of pointsets with pointwise logical operations and
convolutions of their indicator functions, respectively [23]. The later extend
the former by not only providing binary membership tests for set operations,
but also continuous and differentiable membership “grades” that measure
violations of constraints and can be penalized to sensitivity fields for TO
[5, 10].

We overcome the first limitation by defining an analytic generalization of
unsweep to enable incremental co-generation of families of “maximal pairs”
that satisfy the no-collision or (equivalently) containment constraints. We
apply a simple design rule (namely, volume maximization) subject to these
constraints, to encourage persistent contact, i.e., contact maintained
throughout the motion. While this partially overcomes the second limitation
stated above, more sophisticated contact or complementarity measures [7, 8]
can be used in the objective function to enforce persistent contact. It is
important to note that persistent contact depends not only on shape design,
but also on the properties of motion (fixed here) which are not sufficiently
understood. We do not address these issues in this paper to maintain focus on
the co-generation process to incrementally minimize and eliminate collision.
Our approach in this paper is in the same spirit of [23], providing a measure-
theoretic formulation of sweep/unsweep and collision that is differentiable
with respect to the field variables and can thus overcome the third limitation
stated above.

##  3 Collision Measures

In this section, we present a formal modeling framework (Section 3.1) to
define maximal collision-free pairs of solids under given one-parametric
motions, and a discretization scheme (Section 3.2) to make them computable. In
particular, we show that the computations can be factored into offline pre-
computation of a motion-dependent and shape-independent correlation matrix and
online computation of collision measures for specific shapes (i.e., updated
during iterative design) via fast matrix multiplication.

###  3.1 Set-Theoretic Formulation

Let us consider two design domains
$\Omega_{1},\Omega_{2}\subseteq{\mathds{R}}^{\mathrm{d}}$ for the two solid we
aim to co-design ($\mathrm{d}=2,3$). Let us begin from two given initial
designs $S_{1}\subseteq\Omega_{1}$ and $S_{2}\subseteq\Omega_{2}$ representing
the shapes of the two solids at rest (i.e., before applying the motion). Let
$M_{1},M_{2}\in{\mathrm{SE}(\mathrm{d})}$ represent the motions that these
solids experience with respect to a common frame of reference. We restrict our
attention to one-parametric motions defined by the following sets:

$\displaystyle M_{1}\triangleq\big{\\{}\tau_{1}(t)~{}|~{}0\leq\tau\leq
1\big{\\}},$ (3) $\displaystyle
M_{2}\triangleq\big{\\{}\tau_{2}(t)~{}|~{}0\leq\tau\leq 1\big{\\}},$ (4)

where $\tau_{1},\tau_{2}:[0,1]\to{\mathrm{SE}(\mathrm{d})}$ are continuously
time-variant configurations, and can be represented by homogeneous matrices,
vector-quaternion pairs, dual quaternions, etc. The displaced solids at any
given time $t\in[0,1]$ are:

$\displaystyle
S_{1}(t)\triangleq\tau_{1}(t)S_{1}=\big{\\{}\tau_{1}(t){\mathbf{x}}~{}|~{}{\mathbf{x}}\in
S_{1}\big{\\}},$ (5) $\displaystyle
S_{2}(t)\triangleq\tau_{2}(t)S_{2}=\big{\\{}\tau_{2}(t){\mathbf{x}}~{}|~{}{\mathbf{x}}\in
S_{2}\big{\\}}.$ (6)

Without loss of generality, we assume $\tau_{1}(0)=\tau_{2}(0)$ to be
identity, so $S_{1}(0)=S_{1}$ and $S_{2}(0)=S_{2}$, as intended.

To formulate collision measures, it is more convenient to represent the two
pointsets implicitly via binary membership tests, also called indicator or
characteristic functions
${\mathbf{1}}_{S_{1}},{\mathbf{1}}_{S_{2}}:{\mathds{R}}^{3}\to\\{0,1\\}$,
defined generally by:

${\mathbf{1}}_{S}({\mathbf{x}})\triangleq\left\\{\begin{array}[]{ll}1&\text{if}~{}{\mathbf{x}}\in
S,\\\ 0&\text{otherwise}.\end{array}\right.$ (7)

Note that indicator functions are contra-variant with rigid transformations,
i.e., ${\mathbf{1}}_{\tau
S}({\mathbf{x}})={\mathbf{1}}_{S}(\tau^{-1}{\mathbf{x}})$, meaning that a
membership query for a given point against the displaced solid can be computed
by displacing the query point along the inverse trajectory and testing its
membership against the stationary solid.

Let $M=M_{1}^{-1}M_{2}$ stand for the relative motion of $S_{2}$ as observed
from a frame of reference attached to $S_{1}$, noting that by kinematic
inversion, $M^{-1}=M_{2}^{-1}M_{1}$ would represent the relative motion of
$S_{1}$ as observed from a frame of reference attached to $S_{2}$.

Let $M_{1,2}$ (resp. $M_{2,1}$) represent the relative motions of $S_{2}$
(resp. $S_{1}$) as observed from a moving frame of reference attached to
$S_{1}$ (resp. $S_{2}$):

$\displaystyle M_{1,2}\triangleq
M_{1}^{-1}M_{2}=\big{\\{}\tau_{1,2}(t)~{}|~{}0\leq t\leq 1\big{\\}},$ (8)
$\displaystyle M_{2,1}\triangleq
M_{2}^{-1}M_{1}=\big{\\{}\tau_{2,1}(t)~{}|~{}0\leq t\leq 1\big{\\}},$ (9)

where $\tau_{1,2}(t)=\tau_{1}^{-1}(t)\tau_{2}(t)$ and
$\tau_{2,1}(t)=\tau_{2}^{-1}(t)\tau_{1}(t)$. Note also that
$\tau_{1,2}(t)=\tau_{2,1}^{-1}(t)$ hence $M_{1,2}(t)=M_{2,1}^{-1}(t)$.

The displaced solids at any given time $t\in[0,1]$ in the relative frames are:

$\displaystyle
S_{1,2}(t)\triangleq\tau_{2,1}(t)S_{1}=\big{\\{}\tau_{2,1}(t){\mathbf{x}}~{}|~{}{\mathbf{x}}\in
S_{1}\big{\\}},$ (10) $\displaystyle
S_{2,1}(t)\triangleq\tau_{1,2}(t)S_{2}=\big{\\{}\tau_{1,2}(t){\mathbf{x}}~{}|~{}{\mathbf{x}}\in
S_{2}\big{\\}}.$ (11)

To quantify the contribution of every point ${\mathbf{x}}\in{\mathds{R}}^{3}$,
attached to one shape, to its collision with the other shape, we measure the
duration over which the point’s trajectory collides with the latter:

$\displaystyle\mathcal{f}_{2,1}({\mathbf{x}})$
$\displaystyle\triangleq\int_{0}^{1}{\mathbf{1}}_{S_{2,1}(t)}({\mathbf{x}})~{}dt=\int_{0}^{1}{\mathbf{1}}_{S_{2}}(\tau_{2,1}{\mathbf{x}})~{}dt,$
(12) $\displaystyle\mathcal{f}_{1,2}({\mathbf{x}})$
$\displaystyle\triangleq\int_{0}^{1}{\mathbf{1}}_{S_{1,2}(t)}({\mathbf{x}})~{}dt=\int_{0}^{1}{\mathbf{1}}_{S_{1}}(\tau_{1,2}{\mathbf{x}})~{}dt,$
(13)

as illustrated in Fig. 3. To eliminate the contribution of the points that are
outside each shape, we can multiply by the indicator functions of each shape:

$\displaystyle\overline{\mathcal{f}}_{2,1}({\mathbf{x}})$
$\displaystyle\triangleq\int_{0}^{1}{\mathbf{1}}_{S_{2}}(\tau_{2,1}{\mathbf{x}}){\mathbf{1}}_{S_{1}}({\mathbf{x}})~{}dt,$
(14) $\displaystyle\overline{\mathcal{f}}_{1,2}({\mathbf{x}})$
$\displaystyle\triangleq\int_{0}^{1}{\mathbf{1}}_{S_{1}}(\tau_{1,2}{\mathbf{x}}){\mathbf{1}}_{S_{2}}({\mathbf{x}})~{}dt.$
(15)

Figure 3: Local collision measure $\overline{\mathcal{f}}_{2,1}({\mathbf{x}})$
of two query points $\mathbf{x_{1}}$ and
${\mathbf{x}}\in\\{{\mathbf{x}}_{\mathrm{A}},{\mathbf{x}}_{\mathrm{B}}\\}$.
The relative trajectory of each query points is determined by the relative
configurations, $\tau_{1,2}(t)$ of $S_{2}$ to $S_{1}$.

######  Lemma 1.

The two functions in (14) and (15) are implicit representations of the
following unsweeps as their $0-$superlevel sets (up to regularization), i.e.,
the following statement is true “almost everywhere”:333A statement is true
”almost everywhere” for pointsets with respect to a $d-$measure if it be
untrue only over $\mathrm{d}-1$ or lower-dimensional pointsets, which
disappear upon topological regularization [23].

$\displaystyle{\mathbf{x}}\in\mathsf{unsweep}(M_{1,2},\overline{S}_{2})\cap\Omega_{1}~{}\overset{\mathsf{ae}}{\rightleftharpoons}~{}\overline{\mathcal{f}}_{2,1}({\mathbf{x}})=0,$
(16)
$\displaystyle{\mathbf{x}}\in\mathsf{unsweep}(M_{2,1},\overline{S}_{1})\cap\Omega_{2}~{}\overset{\mathsf{ae}}{\rightleftharpoons}~{}\overline{\mathcal{f}}_{1,2}({\mathbf{x}})=0,$
(17)

###### Proof.

The detailed derivation is given in [23]. The idea is that by definition,
points in a given $\mathsf{unsweep}(M,\overline{S})$ are defined by
containment in $\overline{S}$ (i.e., no intersection with $S$) under the
motion $M$. Asserting zero collision measures is a similar criteria, however,
the $\mathrm{d}-$integrals cannot capture the $\mathrm{d}-1$ and lower-
dimensional interference (i.e., contact), hence the equivalence is almost
everywhere. ∎

To derive global measures (a single value for each solid) from the above local
measures, we can integrate them over the respective solids:

$\displaystyle\mathcal{g}_{2,1}$
$\displaystyle\triangleq\int_{S_{1}}\mathcal{f}_{2,1}({\mathbf{x}})~{}d\mu^{\mathrm{d}}[{\mathbf{x}}]=\int_{\Omega_{1}}\overline{\mathcal{f}}_{2,1}({\mathbf{x}})~{}d\mu^{\mathrm{d}}[{\mathbf{x}}],$
(18) $\displaystyle\mathcal{g}_{1,2}$
$\displaystyle\triangleq\int_{S_{2}}\mathcal{f}_{1,2}({\mathbf{x}})~{}d\mu^{\mathrm{d}}[{\mathbf{x}}]=\int_{\Omega_{2}}\overline{\mathcal{f}}_{1,2}({\mathbf{x}})~{}d\mu^{\mathrm{d}}[{\mathbf{x}}],$
(19)

where $d\mu^{\mathrm{d}}[\cdot]$ stands for differential $\mathrm{d}-$measure,
i.e, the area/volume of an infinitesimal 2D/3D region at
${\mathbf{x}}\in{\mathds{R}}^{\mathrm{d}}$. The goal of collision-free co-
design is to find a “maximal pair” of solids $S_{1}\subseteq\Omega_{1}$ and
$S_{2}\subseteq\Omega_{2}$, in a sense that we shall define precisely below,
such that $\mathcal{g}_{S_{2}}=\mathcal{g}_{S_{2}}=0$, which is true iff for
all ${\mathbf{x}}\in{\mathds{R}}^{\mathrm{d}}$,
$\overline{\mathcal{f}}_{2,1}({\mathbf{x}})=\overline{\mathcal{f}}_{1,2}({\mathbf{x}})=0$.

######  Definition 1.

A pair of solids $(S_{1},S_{2})$ are called collision-free under the relative
motion $M_{1,2}=M_{2,1}^{-1}$ if:

$\mu^{\mathrm{d}}[S_{1}\cap
M_{1,2}S_{2}]=0~{}\rightleftharpoons~{}\mu^{\mathrm{d}}[M_{2,1}S_{1}\cap
S_{2}]=0.$ (20)

######  Lemma 2.

A pair of solids $(S_{1},S_{2})$ are collision-free iff the collision measures
are zero, i.e.,
$\mathcal{g}_{2,1}=0~{}\rightleftharpoons~{}\mathcal{g}_{1,2}=0$.

###### Proof.

The above collision criteria implies interference over a
$\mathrm{d}-$dimensional region—while contact over lower ($\mathrm{d}-1$ and
below) dimensional regions is allowed—over a $1-$measurable period of time.
Substituting the local measures in (12) and (13) into the global measures in
(18) and (19) and swapping the space and time integrals yield:

$\displaystyle\mathcal{g}_{2,1}$
$\displaystyle=\int_{0}^{1}\int_{\Omega_{1}}{\mathbf{1}}_{S_{2,1}(t)}({\mathbf{x}}){\mathbf{1}}_{S_{1}}({\mathbf{x}})~{}d\mu^{\mathrm{d}}[{\mathbf{x}}]~{}dt,$
(21) $\displaystyle\mathcal{g}_{1,2}$
$\displaystyle=\int_{0}^{1}\int_{\Omega_{2}}{\mathbf{1}}_{S_{1,2}(t)}({\mathbf{x}}){\mathbf{1}}_{S_{2}}({\mathbf{x}})~{}d\mu^{\mathrm{d}}[{\mathbf{x}}]~{}dt.$
(22)

At a given snapshot $t\in[0,1]$, the inner integrals are nonzero iff
$\mu^{\mathrm{d}}[S_{1}(t)\cap S_{2}(t)]=0$ which, by switching the frame of
reference to the ones attached to the respective solids, is equivalent to
$\mu^{\mathrm{d}}[S_{1}\cap S_{2,1}(t)]=0$ and
$\mu^{\mathrm{d}}[S_{1,2}(t)\cap S_{2}]=0$. The outer integrals are nonzero if
such a nonzero $\mathrm{d}-$measurable interference persists for a
$1-$measurble period of time. ∎

The global collision measures of (21) and (22) are differentiable with respect
to small perturbations of the solids. At a first glance, this might be
counter-intuitive as collisions are understood to be “discrete” events.
However, the measure of collision is continuous with respect to small local
changes in the geometry (e.g., infinitesimal topological inclusions), for a
given relative motion, as discussed in B in greater detail.

######  Definition 2.

A pair of collision-free solids $(S_{1}^{\ast},S_{2}^{\ast})$ under the
relative motion $M_{1,2}=M_{2,1}^{-1}$ is “maximal” if any
$\mathrm{d}-$measurable growth of either solid confined to the design domain
(i.e., $S_{1}^{\ast}\subseteq\Omega_{1}$ and
$S_{2}^{\ast}\subseteq\Omega_{2}$) while keeping the other solid unchanged,
makes them no longer collision-free under the same motion.

More precisely, for every other pair of solids $(S_{1},S_{2})$:

  * •

if $S_{1}^{\ast}\subset S_{1}\subseteq\Omega_{1}$ then $(S_{1},S_{2}^{\ast})$
is not collision-free; and

  * •

if $S_{2}^{\ast}\subset S_{2}\subseteq\Omega_{2}$ then $(S_{1}^{\ast},S_{2})$
is not collision-free.

The maximal pairs can be more formally defined by imposing a partial order
relation $\preccurlyeq$ over the space of all collision-free pairs, where
$(S_{1},S_{2})\preccurlyeq(S_{1}^{\prime},S_{2}^{\prime})$ iff $S_{1}\subseteq
S_{1}^{\prime}$ and $S_{2}=S_{2}^{\prime}$ or $S_{1}=S_{1}^{\prime}$ and
$S_{2}\subseteq S_{2}^{\prime}$. A pair $(S_{1}^{\ast},S_{2}^{\ast})$ is
maximal if
$(S_{1}^{\ast},S_{2}^{\ast})\preccurlyeq(S_{1}^{\prime},S_{2}^{\prime})$
implies $(S_{1}^{\ast},S_{2}^{\ast})=(S_{1}^{\prime},S_{2}^{\prime})$.

Notice that for a maximal pair of collision-free solids
$(S_{1}^{\ast},S_{2}^{\ast})$, if one solid is grown by the slightest bit, the
only way to keep the pair collision-free for the same motion scenario is to
shrink the other solid. This observation provides us with an incremental shape
modification strategy to traverse a family of maximal pairs, akin to walking
along a Pareto front of maximality. If we picture the design space of all
pairs of solids within the domains
$\Omega_{1},\Omega_{2}\subseteq{\mathds{R}}^{\mathrm{d}}$, the maximal front
is a higher-dimensional manifold that bounds the feasible design subspace.

######  Corollary 1.

Given a pair of “initial designs” $S_{1}\subseteq\Omega_{1}$ and
$S_{2}\subseteq\Omega_{2}$ that collide under one-parametric motions
$M_{1},M_{2}\subset{\mathrm{SE}(3)}$, the following two pairs of maximally
collision-free solids can be constructed via unsweep:

$\displaystyle(S_{1},S^{\ast}_{2}),\quad\text{where}~{}S^{\ast}_{2}\triangleq\mathsf{unsweep}(M_{1,2},\overline{S}_{1})\cap\Omega_{2},$
(23)
$\displaystyle(S^{\ast}_{1},S_{2}),\quad\text{where}~{}S^{\ast}_{1}\triangleq\mathsf{unsweep}(M_{2,1},\overline{S}_{2})\cap\Omega_{1},$
(24)

Note that these are extreme cases:

  * •

$S_{2}^{\ast}$ is the maximal solid that does not collide with $S_{1}$,
assuming the shape of $S_{1}$ is fixed.

  * •

$S_{1}^{\ast}$ is the maximal solid that does not collide with $S_{2}$,
assuming the shape of $S_{2}$ is fixed.

However, there are uncountably many other pairs
$(S_{1}^{\prime},S_{2}^{\prime})$ in between the two extreme cases, if we
allow both shapes to change. All of them can be captured by simultaneously
satisfying
$\mathcal{g}_{2^{\prime},1^{\prime}}=\mathcal{g}_{1^{\prime},2^{\prime}}=0$.

We can traverse the family of maximal collision-free pairs
$(S_{1}^{\prime},S_{2}^{\prime})$ using a hyper-parameter $\gamma\in[0,1]$
such that:

$\displaystyle(S_{1}^{\prime},S_{2}^{\prime})_{\gamma}\big{|}_{\gamma:=0}=(S_{1},S^{\ast}_{2}),$
(25)
$\displaystyle(S_{1}^{\prime},S_{2}^{\prime})_{\gamma}\big{|}_{\gamma:=1}=(S^{\ast}_{1},S_{2}),$
(26)

while $(S_{1}^{\prime},S_{2}^{\prime})_{\gamma}$ for $\gamma\in(0,1)$ produces
other maximal pairs in between. Although we cannot provide an explicit formula
for $(S_{1}^{\prime},S_{2}^{\prime})_{\gamma}$, we present an optimization-
based approach to compute the pair in Section 4.

Note that the uncountable family of all maximal pairs cannot be parameterized
with a single (or a number of) parameters. Our goal here is capture a specific
subset of the family with a meaningful hyper-parameter—namely, an enforced
volumetric ratio between the shapes of maximized volume—that serves as a
“knob” to navigate other tradeoffs for design space exploration.

###  3.2 Computational Discretization

The above formulation is representation-agnostic and can be used with a
variety of different representation schemes (e.g., B-reps, mesh, and voxels),
as long as they support computation of $\mathrm{d}-$integrals—and if they do
not, one can approximate $\mathrm{d}-$integrals by sampling the solids using
quadrature rules, as long as the representation schemes support basic
membership classification queries. Here, we select a simple asymmetric
discretization approach to illustrate the practicality of our formulation.

To enable finite representation of solids and digital computation of collision
measures, let us approximate the indicator function of a stationary solid
$S\subseteq{\mathds{R}}^{\mathrm{d}}$ using a finite volume (FV) scheme over a
grid, meaning that each grid cell is associated with densities (i.e.,
$\mathrm{d}-$measure fractions) of the solid entrapped within the cell,
denoted by $\rho_{S,i}\in[0,1]$ for $i=1,2,\ldots,n$. For a moving solid
$S\subseteq{\mathds{R}}^{\mathrm{d}}$, on the other hand, let us use a finite
sample (FS) scheme over the same grid, meaning that each cell center is
associated with the same density values like a “lumped” measure. A rigorous
application of these discretization schemes to approximate the collision
integrals by sums in a $\mathrm{d}-$measure-preserving fashion over any subset
of the grid cell is presented in A. The resulting sums are:

$\displaystyle\mathcal{g}_{2,1}$
$\displaystyle\approx\sum_{i_{1}=1}^{n_{1}}\sum_{i_{2}=1}^{n_{2}}\rho_{S_{1},i_{1}}\rho_{S_{2},i_{2}}w^{1,2}_{i_{1},i_{2}},$
(27) $\displaystyle\mathcal{g}_{1,2}$
$\displaystyle\approx\sum_{i_{1}=1}^{n_{1}}\sum_{i_{2}=1}^{n_{2}}\rho_{S_{1},i_{1}}\rho_{S_{2},i_{2}}w^{2,1}_{i_{2},i_{1}},$
(28)

where the weights $w^{1,2}_{i_{1},i_{2}}$ and $w^{2,1}_{i_{2},i_{1}}$ mean the
following:

  * •

$w^{1,2}_{i_{1},i_{2}}$ measures how long the grid vertex
${\mathbf{x}}_{i_{2}}\in\Omega_{2}$ of the moving grid (attached to
$S_{2}\subseteq\Omega_{2}$) stays within the $\mathrm{d}-$cell $C_{i_{1}}$
under the motion $\tau_{1,2}(t)$ for $t\in[0,1]$.

  * •

$w^{2,1}_{i_{2},i_{1}}$ measures how long the grid vertex
${\mathbf{x}}_{i_{1}}\in\Omega_{1}$ of the moving grid (attached to
$S_{1}\subseteq\Omega_{1}$) stays within the $\mathrm{d}-$cell $C_{i_{2}}$
under the motion $\tau_{2,1}(t)$ for $t\in[0,1]$.

For further details on the FV and FS discretization schemes, see A.

The finite approximations of the collision measures in (27) and (28) can be
written as matrix equations:

$\displaystyle\mathcal{g}_{2,1}$
$\displaystyle\approx\big{[}\rho_{S_{1},i_{1}}\big{]}^{\mathrm{T}}\big{[}w^{1,2}_{i_{1},i_{2}}\big{]}\big{[}\rho_{S_{2},i_{2}}\big{]},$
(29) $\displaystyle\mathcal{g}_{1,2}$
$\displaystyle\approx\big{[}\rho_{S_{2},i_{2}}\big{]}^{\mathrm{T}}\big{[}w^{2,1}_{i_{2},i_{1}}\big{]}\big{[}\rho_{S_{1},i_{1}}\big{]},$
(30)

The two arrays $[\rho_{S_{1},i_{1}}]_{n_{1}\times 1}$ and
$[\rho_{S_{2},i_{2}}]_{n_{2}\times 1}$ are discrete representations of the two
solids, i.e., the design variables. Importantly, the two matrices
$[w^{1,2}_{i_{1},i_{2}}]_{n_{1}\times n_{2}}$ and
$[w^{2,1}_{i_{2},i_{1}}]_{n_{2}\times n_{1}}$ do not depend on the designs.
The weights can be viewed as pairwise correlations between primal grid nodes
of a moving grid and dual grid cells of a stationary grid, which depend solely
on the relative motion of the grids and the grid structure. The matrices can
thus be pre-computed offline and reused in iterative design of the two solids.
This property is critical for computational tractability of iterative design
optimization in Section 4, given that computing the collision measures for
arbitrarily complex shapes and motions can be prohibitive in an iterative loop
running hundreds or thousands of iterations.

Figure 4 illustrates how the stationary and moving solids are discretized via
FV and FS, respectively, and how the correlation matrix entries are computed
in practice.

Figure 4: The stationary solid (in this case, $S_{1}$) is discretized via FV
while the moving solid (in this case, $S_{2}$) is discretized via FS. The
correlation matrix entry $w^{2,1}_{i_{2},i_{1}}$ is computed by counting the
number of times the discretized trajectory of the cell center at $i_{1}$ of
the stationary grid stays within the cell at $i_{2}$ of the moving grid.

To enable co-design of $[\rho_{S_{1},i_{1}}]_{n_{1}\times 1}$ and
$[\rho_{S_{2},i_{2}}]_{n_{2}\times 1}$ subject to collision avoidance
constraints, the violation of such constraints (i.e., the collision measures)
must be differentiated with respect to the design variables. The resulting
discrete sensitivity fields associated with the respective grids are computed
using a chain rule:

$\displaystyle\left[\frac{\partial\mathcal{g}_{2,1}}{\partial\rho_{S_{1},i_{1}}}\right]$
$\displaystyle\approx\big{[}w^{1,2}_{i_{1},i_{2}}\big{]}\big{[}\rho_{S_{2},i_{2}}\big{]},~{}\left[\frac{\partial\mathcal{g}_{2,1}}{\partial\rho_{S_{2},i_{2}}}\right]\approx\big{[}w^{1,2}_{i_{1},i_{2}}\big{]}^{\mathrm{T}}\big{[}\rho_{S_{1},i_{1}}\big{]},$
(31)
$\displaystyle\left[\frac{\partial\mathcal{g}_{1,2}}{\partial\rho_{S_{1},i_{1}}}\right]$
$\displaystyle\approx\big{[}w^{2,1}_{i_{2},i_{1}}\big{]}\big{[}\rho_{S_{1},i_{1}}\big{]},~{}\left[\frac{\partial\mathcal{g}_{1,2}}{\partial\rho_{S_{2},i_{2}}}\right]\approx\big{[}w^{2,1}_{i_{2},i_{1}}\big{]}^{\mathrm{T}}\big{[}\rho_{S_{2},i_{2}}\big{]},$
(32)

See B for a definition of the topological sensitivity fields (TSFs), defined
for the continuum (i.e., in terms of the indicator functions), corresponding
to $\epsilon\to 0^{+}$.

##  4 Generation of collision-free Geometries

In this section, we present an optimization based formulation (Section 4.1)
that allows both collision-free solids to evolve in a controlled manner, using
the collision measures developed in the previous section.

###  4.1 Optimization Problem Formulation

Using collision measures as constraints is not sufficient to define a well-
posed optimization problem, because every subset of the maximal solids is also
collision-free (despite not maximal). In practice, attempting to minimize the
collision measures by gradient-descent optimization can lead to useless
designs where the solids are excessively shrunk to avoid collision without
maintaining contact. Figure 5 shows an example in 2D; we begin with a
stationary disk and a rotating disk around its center, and can satisfy
collision avoidance with uncountably many solutions, adding/removing material
from either disk. However, only two of them (I and II) are maximal pairs while
(III) is not. Therefore, additional constraints are needed to achieve
maximality.

Figure 5: Two initial designs $S_{1}$ and $S_{2}$ are shown (top left), where
$S_{2}$ moves relative to $S_{1}$ by a rotation about the center of $S_{2}$.
Three collision-free pairs of solids are shown for this motion, where in (a)
and (b) the pairs are maximal while in (c) they are not. Note that $S_{2}$ is
the entire disk on the right including the white and teal regions.

First, observe that we can partition the initial designs
$S_{1}\subseteq\Omega_{1}$ and $S_{2}\subseteq\Omega_{2}$ into two subdomains
each; namely,

  * •

The initially colliding subdomains $\widehat{S}_{1}$ and $\widehat{S}_{2}$
are:

$\displaystyle\widehat{S}_{1}\triangleq\big{\\{}{\mathbf{x}}\in
S_{1}~{}|~{}\overline{\mathcal{f}}_{2,1}({\mathbf{x}})>0\big{\\}},$ (33)
$\displaystyle\widehat{S}_{2}\triangleq\big{\\{}{\mathbf{x}}\in
S_{2}~{}|~{}\overline{\mathcal{f}}_{1,2}({\mathbf{x}})>0\big{\\}},$ (34)

  * •

The initially collision-free subdomains $\widetilde{S}_{1}$ and
$\widetilde{S}_{2}$ are:

$\displaystyle\widetilde{S}_{1}\triangleq\big{\\{}{\mathbf{x}}\in
S_{1}~{}|~{}\overline{\mathcal{f}}_{2,1}({\mathbf{x}})=0\big{\\}},$ (35)
$\displaystyle\widetilde{S}_{2}\triangleq\big{\\{}{\mathbf{x}}\in
S_{2}~{}|~{}\overline{\mathcal{f}}_{1,2}({\mathbf{x}})=0\big{\\}}.$ (36)

Note that $S_{1}=\widehat{S}_{1}\cup\widetilde{S}_{1}$ and
$S_{2}=\widehat{S}_{2}\cup\widetilde{S}_{2}$. We only need to eliminate the
collisions by modifying the latter.

The simplest possible way to make the optimization problem well-posed is to
pick an objective function that pushes the collision avoidance constraints to
become active, i.e., the resulting pairs of solids become maximally collision-
free. We observed that maximizing the $\mathrm{d}-$measure of both design
subdomains (e.g., area in 2D and volume in 3D) is effective in many
scenarios—although more sophisticated objective functions such as contact [7]
and complementarity [8] measures can be adopted in future work. Moreover, to
control the choice of maximal pairs from a one-parametric family using a
meaningful hyper-parameter $\gamma\in[0,1]$, we add another constraint that
enforces a $\mathrm{d}-$measure ratio between the two solids (37d). The
optimization problem is thus formulated as finding
$\widehat{S}_{1}\subseteq\Omega_{1}$ and $\widehat{S_{2}}\subseteq\Omega_{2}$
to:

$\displaystyle\mathop{\text{Maximize}}\limits_{\widehat{S}_{1},\widehat{S_{2}}}\quad$
$\displaystyle
f(\widehat{S}_{1},\widehat{S}_{2})=\mu^{\mathrm{d}}[\widehat{S}_{1}]+\mu^{\mathrm{d}}[\widehat{S}_{2}],$
(37a) subject to $\displaystyle{\mathcal{g}}_{\widehat{2},\widehat{1}}=0,$
(37b) $\displaystyle{\mathcal{g}}_{\widehat{1},\widehat{2}}=0,$ (37c)
$\displaystyle
g(\widehat{S}_{1},\widehat{S}_{2})=\gamma\mu^{\mathrm{d}}[\widehat{S}_{1}]-(1-\gamma)\mu^{\mathrm{d}}[\widehat{S}_{2}]=0$
(37d)

For gradient-descent optimization, we define a Lagrangian from the general
formulation:

$\mathcal{L}(\widehat{S}_{1},\widehat{S}_{2})=f(\widehat{S}_{1},\widehat{S}_{2})+\lambda_{1}\mathcal{g}_{\widehat{2},\widehat{1}}+\lambda_{2}\mathcal{g}_{\widehat{1},\widehat{2}}+\lambda_{3}g(\widehat{S}_{1},\widehat{S}_{2}).$
(38)

Using the discretization scheme presented in Section 3.2, both
$\widehat{S}_{1}\subseteq\Omega_{1}$ and $\widehat{S_{2}}\subseteq\Omega_{2}$
can be represented by the density arrays
$[\rho_{\widehat{S}_{1}}]_{n_{1}\times 1}$ and
$[\rho_{\widehat{S}_{2}}]_{n_{2}\times 1}$, with pairwise collisions,
$\big{[}\widehat{w}^{1,2}_{i_{1},i_{2}}\big{]}$ and
$\big{[}\widehat{w}^{2,1}_{i_{2},i_{1}}\big{]}$ hence:

$\displaystyle f(\widehat{S}_{1},\widehat{S}_{2})$
$\displaystyle\approx\epsilon^{\mathrm{d}}\big{\|}\rho_{\widehat{S}_{1}}\big{\|}_{1}+\epsilon^{\mathrm{d}}\big{\|}\rho_{\widehat{S}_{2}}\big{\|}_{1},$
(39) $\displaystyle g(\widehat{S}_{1},\widehat{S}_{2})$
$\displaystyle\approx\gamma\epsilon^{\mathrm{d}}\big{\|}\rho_{\widehat{S}_{1}}\big{\|}_{1}-(1-\gamma)\epsilon^{\mathrm{d}}\big{\|}\rho_{\widehat{S}_{2}}\big{\|}_{1},$
(40)

where $\|\cdot\|_{1}$ is the $L_{1}-$norm, i.e.,
$\|\rho_{\widehat{S}_{1}}\|_{1}$ and $\|\rho_{\widehat{S}_{2}}\|_{1}$ are sum
of non-negative density values in each array, approximating the
$\mathrm{d}-$measure of the solids they represent:

$\displaystyle\big{\|}\rho_{\widehat{S}_{1}}\big{\|}_{1}$
$\displaystyle\triangleq\sum_{i_{1}=1}^{n_{1}}\big{|}\rho_{\widehat{S}_{1},i_{1}}\big{|}=\sum_{i_{1}=1}^{n_{1}}\rho_{\widehat{S}_{1},i_{1}},$
(41) $\displaystyle\big{\|}\rho_{\widehat{S}_{2}}\big{\|}_{1}$
$\displaystyle\triangleq\sum_{i_{2}=1}^{n_{2}}\big{|}\rho_{\widehat{S}_{2},i_{2}}\big{|}=\sum_{i_{2}=1}^{n_{2}}\rho_{\widehat{S}_{2},i_{2}}.$
(42)

Without loss of generality, we can eliminate $\epsilon:=1$ from all equations,
given that our approach is scale-agnostic. Substituting for the collision
measures in (29) and (30) as well as the above norms in (38), we obtain:

$\displaystyle\mathcal{L}(\widehat{S}_{1},\widehat{S}_{2})$
$\displaystyle\approx(1+\lambda_{3}\gamma)\big{\|}\rho_{\widehat{S}_{1}}\big{\|}_{1}+(1-\lambda_{3}(1-\gamma))\big{\|}\rho_{\widehat{S}_{2}}\big{\|}_{1}$
$\displaystyle+\lambda_{1}\big{[}\rho_{\widehat{S}_{1},i_{1}}\big{]}^{\mathrm{T}}\big{[}\widehat{w}^{1,2}_{i_{1},i_{2}}\big{]}\big{[}\rho_{\widehat{S}_{2},i_{2}}\big{]}$
$\displaystyle+\lambda_{2}\big{[}\rho_{\widehat{S}_{2},i_{2}}\big{]}^{\mathrm{T}}\big{[}\widehat{w}^{2,1}_{i_{2},i_{1}}\big{]}\big{[}\rho_{\widehat{S}_{1},i_{1}}\big{]}.$
(43)

All of the quantities in the above Lagrangian and its partial derivatives with
respect to design variables can be readily computed through linear-algebraic
operations, which scale well for parallel computing on CPU/GPU, although we
have not implemented such parallelization in this paper.

###  4.2 Iterative Optimization Algorthm

To solve the optimization problem an iterative approach is utilized where the
method of moving asymptotes (MMA) [29] is utilized. Algorithm 1 provides a
description of the approach.

Algorithm 1 Co-generation of maximal collision-free solids

procedure Co-generation($S_{1},S_{2},M_{1,2},M_{2,1},\gamma$)

Initialize $[\rho_{S_{1}}]\leftarrow[{\mathbf{1}}_{S_{1}}]$

Initialize $[\rho_{S_{2}}]\leftarrow[{\mathbf{1}}_{S_{2}}]$

Initialize $\Delta\leftarrow\infty$

Initialize $iter\leftarrow 0$

$(\big{[}w^{1,2}_{i_{1},i_{2}}\big{]},\big{[}w^{2,1}_{i_{2},i_{1}}\big{]})\leftarrow\textsc{PrecompCo}(S_{1},S_{2},M_{1,2},M_{2,1})$

$([\widehat{\rho}],\big{[}\widehat{w}^{1,2}_{i_{1},i_{2}}\big{]},\big{[}\widehat{w}^{2,1}_{i_{2},i_{1}}\big{]})\leftarrow\textsc{ColVox}([{\rho}],\big{[}w^{1,2}_{i_{1},i_{2}}\big{]},\big{[}w^{2,1}_{i_{2},i_{1}}\big{]})$

while $\Delta>\delta$ and $iter<l$ do

$f\leftarrow\textsc{\text{Evaluate}}([\widehat{\rho}])$ $\triangleright$ Obj.
func.

$\frac{\partial\mathcal{L}_{[\widehat{\rho}]}}{\partial\widehat{\rho}}\leftarrow\textsc{\text{Gradient}}([f,\widehat{\rho},\big{[}\widehat{w}^{1,2}_{i_{1},i_{2}}\big{]},\big{[}\widehat{w}^{2,1}_{i_{2},i_{1}}\big{]})$$\triangleright$
Sens.

$[\widehat{\rho}^{\text{new}}]\leftarrow\textsc{\text{Update}}([\widehat{\rho}],\frac{\partial\mathcal{L}_{[\widehat{\rho}]}}{\partial\widehat{\rho}})$
$\triangleright$ Gradient Update

$\Delta\leftarrow\textsc{Integrate}([\widehat{\rho}^{\text{new}}]-[\widehat{\rho}])$
$\triangleright$ Vol. diff.

$iter\leftarrow iter+1$ $\triangleright$ Iter. counter

$[\widehat{\rho}]\leftarrow[\widehat{\rho}^{\text{new}}]$ $\triangleright$ For
next iteration

end while

return $[\widehat{\rho}]$

end procedure

With the optimization formulated, we provide various examples of co-generation
of collision-free solids.

##  5 Results & Discussion

To demonstrate the efficacy of the approach, we provide a variety of 2D and 3D
examples. All examples are run on a desktop machine with Intel® CoreTM
i7-9800X CPU with 16 processors running at 3.8 GHz, 32 GB of host memory, and
an NVIDIA® GeForce® GTX 1080 GPU with 2,560 CUDA cores and 8 GB of device
memory.

###  5.1 2D Counter Rotating Squares

The first example for the algorithm is two counter rotating squares where the
angular position of the squares are provided by $\theta_{1}\in[0,+2\pi)$ and
$\theta_{2}\in[0,-2\pi)$ (counter-clockwise) where the motion specified by
$\theta_{1}=-\theta_{2}=\omega t$ with constant angular velocity $\omega>0$.
The dimensions and initial angular positions of the squares are shown in Fig.
6 (top). Additionally, the normalized local collision measures are provided in
Fig. 6 (bottom) to show how initially the two shapes collide. The temporal
resolution is $500$ time steps and the spatial resolution is $400\times 400$
pixels.

Figure 6: Two initial design (top), $S_{1}$ and $S_{2}$, where each solid
rotate around their centers with the same angular speed in opposite
directions. The local collision measure for each domain is shown (bottom).

The co-generation procedure was applied for $\gamma$ varying from 0 to 1 in
0.1 increments. To summarize the results, Fig. 7 provides a plot of the sum of
the $\mathrm{d}-$measure of each co-generated solid as a function of $\gamma$.
Selected solutions co-generated by the procedure are included as well. The
solids are thresholded at $[{\rho}_{S_{1},i_{1}}],[{\rho}_{S_{2},i_{2}}]>0.5$.

Figure 7: Sum of the $\mathrm{d}-$measures of the two collision-free solids as
a function of hyper-parameter $\gamma$.

Convergence plots are shown in Fig. 8 for $\gamma=0.4$ and $\gamma=0.5$ to
demonstrate the convergence of the measure during incremental optimization.

Figure 8: Convergence of summed $\mathrm{d}-$measures of the two collision-
free solids for $\gamma=0.4$ (left) and $\gamma=0.5$ (right).

Due to the symmetry of the trajectories and domains, the total weighted area
as a function of $\gamma$ is symmetric about $\gamma=0.5$. The solids also
reflect this symmetry. Additionally, the results of $\gamma=0$ and $\gamma=1$
precisely match the results that would be obtained if the unsweep operations
in (23) and (24) were performed.

###  5.2 2D Cam and Follower

The second example is a cam/follower system. The cam, which is initially a
square of length $L$, is prescribed to rotate $2\pi$ radians about its center
(the origin $O$) at a constant angular velocity, while the follower translates
in the vertical direction. The vertical position $y_{\mathrm{F}}$ of the
center of the follower is prescribed as a function of the angular position
$\theta_{C}$ of the cam:

$y_{\mathrm{F}}={{3L}\over{4}}+{{L}\over{8}}{\cos(2\theta_{\mathrm{C}})}.$
(44)

The dimensions and initial positions are shown in Fig. 9. Additionally in Fig.
9 the normalized local collision measures are provided to show how the two
shapes collide before incrementally changing their shapes. The plots of local
collision measures are offset and rescaled for clarity and do not correspond
to the initial physical locations.

Figure 9: Two initial designs, cam ($S_{\mathrm{C}}$) and follower
($S_{\mathrm{F}}$) (left). The local collision measure for each domain is
shown (right).

The temporal resolution is 1,000 time steps while the spatial resolutions are
$400\times 400$ for the cam and $400\times 200$ for the follower. This problem
is particularly interesting because the unsweep operations result in a loss of
contact at certain times in the trajectory as shown in Fig. 10. Our co-
generation process produces persistent contact to ensure functionality,
detailed below.

Figure 10: Illustration of how the unsweep operation can result in a loss of
contact for the cam and follower.

The co-generation procedure was applied for $\gamma$ varying from 0 to 1 in
0.1 increments. Fig. 11 provides a plot of the average minimum distance
between the two solutions as they move through their trajectories as a
function of $\gamma$. As discussed, certain $\gamma$ values result in a loss
of contact and therefore would not be functional. But in the range of
approximately $\gamma\in(0.7,0.9)$ the solutions maintains contact even though
it is not an explicit constraint. Remember that our approach does not
guarantee persistent contact, as (a) it may not even be possible for any
maximal pair for the given motion; (b) the proper shapes may not be obtainable
by the specific subset of maximal pairs parameterized by $\gamma$ and computed
by our incremental procedure. In particular, maximality in terms of volumes
and collision avoidance may imply contact in one configuration, but not
persistent contact throughout the whole motion, as shown in the results of
Fig. 11. Nevertheless, our observation is that co-generation is a necessary
approach to create functional parts. Future work should focus on studying
shape and motion properties with regards to persistent contact.

Figure 11: Average minimum distance between cam and follower as a function of
$\gamma$ which illustrates co-generation is necessary to maintain contact
during motion.

To demonstrate the improved contact via co-generation, the solids generated
when $\gamma=0.8$ are shown in Fig. 12 at various positions during its
trajectory. The cam and follower never lose contact during motion which is a
necessary condition for functionality.

Figure 12: Cam and follower systems co-generated with $\gamma=0.8$ shown at
various configurations during the motion.

The discretized sensitivity fields defined in (31) and (32) are plotted in
Fig. 13 during the first iteration of the optimization loop.

Figure 13: Plots showing the four discrete sensitivity fields for the initial
2D cam and follower example.

###  5.3 3D Cam and Follower

The third example is a 3D cam/follower system. The cam rotates around the
$z-$axis by an angle $\theta_{\mathrm{C}}$ at a constant angular velocity.
Simultaneously, the follower rotates around the $x-$axis with by an angle
$\theta_{\mathrm{F}}$ where
$\theta_{\mathrm{F}}=\nicefrac{{1}}{{2}}|\sin(\theta_{\mathrm{C}})|$. The
rotation matrices describing the motion of the cam and follower,
$R_{\mathrm{C}}$ and $R_{\mathrm{C}}$, respectively, are provided below:

$\displaystyle R_{\mathrm{C}}$
$\displaystyle=\begin{bmatrix}~{}\cos(\theta_{\mathrm{C}})&~{}\sin(\theta_{\mathrm{C}})&0\\\
-\sin(\theta_{\mathrm{C}})&~{}\cos(\theta_{\mathrm{C}})&0\\\ 0&0&1\\\
\end{bmatrix},$ (45) $\displaystyle R_{\mathrm{F}}$
$\displaystyle=\begin{bmatrix}1&0&0\\\
0&~{}\cos(\theta_{\mathrm{F}})&~{}\sin(\theta_{\mathrm{F}})\\\
0&-\sin(\theta_{\mathrm{F}})&~{}\cos(\theta_{\mathrm{F}})\\\ \end{bmatrix},$
(46)

The design domain dimensions at initial configurations are shown in Fig. 14.
The temporal resolution is 5,000 time steps, while the spatial resolutions are
$105\times 105\times 105$ for the cam and $105\times 105\times 140$ for the
follower. Each domain has over a million voxels. Similar to the 2D
cam/follower system, if the unsweep operation was used to generate the
collision-free shapes, a loss of contact would occur. By co-generating the
solids, on the other hand, contact can be maintained as shown in Fig. 14,
where the cam and follower are rendered in various relative configurations
during the motion for $\gamma=0.5$. In addition, a convergence plot is
provided for $\gamma=0.5$ in Fig. 15.

Figure 14: Initial solids for a spherical cam and follower. Figure 15:
Convergence of summed $\mathrm{d}-$measures of the two collision-free solids
for $\gamma=0.5$.

###  5.4 3D Bolt and Nut

The final example is a 3D bolt and nut pair. The dimensions and initial
positions are shown in Fig. 16, where $\phi_{\mathrm{S}}$ is the rotation
angle of the bolt around the $z-$axis with a prescribed pitch of
$P_{\mathrm{S}}=\nicefrac{{L}}{{5}}$, while the nut is stationary. The
homogeneous transformation matrix $T_{\mathrm{S}}$, describing the relative
screw motion is:

$R_{\mathrm{S}}=\begin{bmatrix}~{}\cos(\phi_{\mathrm{S}})&~{}\sin(\phi_{\mathrm{S}})&0&0\\\
-\sin(\phi_{\mathrm{S}})&~{}\cos(\phi_{\mathrm{S}})&0&0\\\
0&0&1&{{-L}\over{10\pi}}\phi\\\ 0&0&0&1\\\ \end{bmatrix},\\\ $ (47)

where the translation ${{-L}\over{10\pi}}\phi_{\mathrm{S}}$ along and rotation
$\phi_{\mathrm{S}}\in[0,8\pi)$ around the $z-$axis are linearly related.
During the motion, the bolt makes 4 full turns and moves by
$4P_{\mathrm{S}}=0.8L$. The temporal resolution is 5,000 time steps, while the
spatial resolutions are $50\times 50\times 150$ for the bolt and $100\times
100\times 50$ for the nut. One particular solution in which $\gamma=0.2$ is
provided in Fig. 16.

Figure 16: Inital design domains (left) and collision-free bolt and nut
generated (right) with threads naturally emerging from co-design.

The automatic co-generation of threads and convergence to a bolt and nut pair,
starting from two colliding cubic regions, is a significant result. The
threads appear naturally and automatically to generate maximal collision-free
pairs for a given screw motion, leading to a bolt and nut pair with the same
pitch as that of the given screw motion.

Table 1: Summary of computation time for various examples provided in Section 5. Case Study | Spatial Res. (1) | Spatial Res. (2) | Temporal Res. | Mean Comp. Time (sec) | Max Iters.
---|---|---|---|---|---
2D Gear/Gear | 400$\times$400 | 400$\times$400 | 500 | 705.2 | 150
2D Cam/Follower | 400$\times$400 | 400$\times$400 | 1,000 | 949.1 | 200
3D Cam/Follower | 105$\times$105$\times$105 | 105$\times$105$\times$140 | 5,000 | 7,890 | 500
3D Bolt/Nut | 50$\times$50$\times$50 | 100$\times$100$\times$50 | 5,000 | 4,201 | 350

However, note that the bolt and nut pair is not the only locally optimal
solutions for a screw motion according to the formulation in (37). A
cylindrical peg and hole pair, for example, would satisfy all the enforced
constraints, including maximal collision-avoidance, maximal (space-filling)
volume, and desired volume ratio, while ensuring persistent contact as a side-
effect. However, there are infinitely many screw and nut pairs (of different
thread profiles) that could satisfy the constraints with the same maximal
volume. Hence, it appears that finding one of these threaded solids is more
likely than finding the degenerate extreme case of cylindrical peg and hole
(i.e., no threads). It appears that our approach converged to one such local
extremum of the optimization problem.

It is also worth noting that this procedure certainly does not generate all
possible maximally collision-free pairs. Additional constraints or an entirely
different set of objective functions and constraints—e.g., maximizing contact
[7] or complementarity [8] measures as opposed to total volume and volume
ratio—may be used to identify more maximal pairs with better form, fit, or
function.

Table 1 summarizes the presented results. The spatial resolution of both
design domains, the temporal resolution of the motions, the mean computation
time for each optimization run, and the maximum number of optimization loop
iterations are presented. Our numerical procedure required running times on
the order of minutes and hours, indicating the scalability of the approach for
high-resolution spatio-temporal discretizations.

###  5.5 Summary of Results & Discussion

By providing various examples, we demonstrated the ability of our approach to
generate pairs of collision-free solids in both 2D and 3D, discretized by
immersing them into a pair of Cartesian grids using FV and FS schemes. By
precomputing the sparse pairwise correlation matrix for a pair of grids in
relative motion, we demonstrated the scalability of the approach. The most
computationally demanding task in the design loop is cast into a pre- and
post-multiplication of the correlation matrix with the design variables
(namely, the density arrays for each solid). The results indicate that co-
generation is essential, not only for generating a broad family of collision-
free solids to enable downstream design flexibility, but also to find well-
fitting solids that maintain contact throughout the motion.

There are a few seemingly arbitrary choices that may affect the outcome of our
approach. While the choice of initial solids can impact the final design in
the extreme cases (e.g., when they are too small to collide resulting in an
ill-posed optimization problem) the initial shapes appear less consequential
than the motion and choice of the hyper-parameter (i.e., volume ratio)
$\gamma$. There are uncountable pairs of continuum solids and exponentially
many pairs of corresponding discrete (e.g., voxelized) solids that are
maximal, among which we have chosen a one-dimensional subset parameterized by
$\gamma$. A natural extension of the current work is to explore higher-
dimensional subsets with more parameters, providing multiple “knobs” to
designers to adjust additional/implicit design requirements while satisfying
maximal collision-free contact.

Hopefully, our work will inspire formulating better objective functions and
constraints that encode complex requirements such as physics, contact, and
complementarity. A particularly promising direction is to combine the maximal
collision avoidance with physics-based performance criteria (e.g., computed by
FEA for a given set of loading conditions on the parts), manufacturability
[10, 12], and other constraints [30] using the integration framework in [5].

##  6 Conclusion

In this paper we introduced an automated and scalable procedure for generative
co-design of pairs of solids under one-parametric motions to satisfy collision
avoidance or (equivalently) containment constraints. Our automated procedure
designs a family of “maximal pairs” of collision-free solids simultaneously.
Our procedure can design shapes that maintain contact (when possible) while
avoiding collisions, using a hyper-parameter to control the relative volume of
their modified regions. Unlike existing methods that require making arbitrary
restrictions such as parameterizing the shapes for common/simple motions or
fixing one shape and computing the other via unsweep, our design process
avoids premature decisions and unnecessary restrictions that lead to
suboptimal solutions. We also show that unsweep of each solid against the
other is subsumed as a special case of our operation at the two extremes of
the hyper-parameter spectrum.

While collision avoidance is not sufficient to ensure persistent contact, our
results suggest that incremental co-design to reduce collision while
maximizing volume (or perhaps other objective functions) is a viable route for
designing mechanisms in which contact and complementarity are critical for
function, such as cams/followers and gear trains. More effective contact or
complementarity measures can likely be employed, using similar measure-
theoretic approaches used to quantify collision, to replace the simple volume
maximization. Additionally, the formulation of the collision measures as
differentiable fields enables seamless integration into gradient-based design
optimization approaches such as TO. In the context of design for assembly,
this could broaden the scope of TO and generative design from part-level to
assembly-level design, in which parts and sub-assemblies are simultaneously
evolved for physics-based performance, manufacturing, and other objectives
while ensuring parts are collision-free.

It should be noted that when co-designing solids within larger assemblies
(e.g., 100s of parts), one has to deal with combinatorial explosion of
collision analysis among pairs, triplets, and so on. Future work in this
direction should explore efficient ways to avoid such exhaustive computations.
In the worst case scenario, the initial design domain of each component can be
selected largely enough for all parts in an assembly to collide with each
other. Computing all of the collision matrices on top of simultaneously
optimizing all parts would be computationally intractable. One possible
mitigation is to come up with domain restriction strategies or spatial
sorting/hashing to determine components that are collision-free without
computing all motions. In addition to prohibitive running times, precomputing
the correlation matrices will be memory-intensive. While the matrices are
often sparse, it may not be the case for exceptionally high spatial
resolutions and complex motions leading to large numbers of point-voxel
collisions in each time-step. Finally the use of symmetries based on the
motions should be explored to reduce computational time and memory costs.

## Acknowledgments

The authors are thankful to Amir Mirzendehdel for his insights on topology
optimization and Saigopal Nelaturi for his support of this project under the
Design for Assembly research program at PARC.

## References

## References

  * [1] G. Boothroyd, Design for manufacture and assembly: The Boothroyd-Dewhurst experience, in: Design for X, Springer, 1996, pp. 19–40.
  * [2] L. Chong, S. Ramakrishna, S. Singh, A review of digital manufacturing-based hybrid additive manufacturing processes, The International Journal of Advanced Manufacturing Technology 95 (5) (2018) 2281–2300.
  * [3] J. Sokolowski, J. P. Zolésio, Introduction to shape optimization, Springer, 1992, pp. 5–12.
  * [4] O. Sigmund, K. Maute, Topology optimization approaches, Structural and Multidisciplinary Optimization 48 (6) (2013) 1031–1055.
  * [5] A. M. Mirzendehdel, M. Behandish, S. Nelaturi, Exploring feasible design spaces for heterogeneous constraints, Computer-Aided Design 115 (2019) 323–347.
  * [6] M. Lysenko, Fourier collision detection, The International Journal of Robotics Research 32 (4) (2013) 483–503.
  * [7] M. Lysenko, V. Shapiro, Effective contact measures, Computer-Aided Design 70 (2016) 134–143.
  * [8] M. Behandish, H. T. Ilieş, Shape complementarity analysis for objects of arbitrary shape, arXiv preprint arXiv:1712.00238.
  * [9] R. L. Norton, Design of Machinery: An Introduction to the Synthesis and Analysis of Mechanisms and Machines, McGraw-Hill Higher Education, 2008.
  * [10] A. M. Mirzendehdel, M. Behandish, S. Nelaturi, Topology pptimization with accessibility constraint for multi-axis machining, Computer-Aided Design 122 (2020) 102825.
  * [11] A. M. Mirzendehdel, M. Behandish, S. Nelaturi, Optimizing build orientation for support removal using multi-axis machining, Computers & Graphics 99 (2021) 247–258.
  * [12] A. M. Mirzendehdel, M. Behandish, S. Nelaturi, Topology optimization for manufacturing with accessible support structures, Computer-Aided Design 142 (2022) 103117.
  * [13] D. Chakravarthi, Two cylinder engine-mechanism, https://grabcad.com/library/two-cylinder-engine-mechanism-1 (2022).
  * [14] H. T. Ilieş, V. Shapiro, On shaping with motion, J. Mech. Des. 122 (4) (2000) 567–574.
  * [15] H. T. Ilieş, V. Shapiro, A class of forms from function: The case of parts moving in contact, Research in Engineering Design 13 (3) (2002) 157–166.
  * [16] H. T. Ilieş, V. Shapiro, The dual of sweep, Computer-Aided Design 31 (3) (1999) 185–201.
  * [17] H. T. Ilieş, V. Shapiro, UNSWEEP: Formulation and computational properties, in: Proceedings of the 4th ACM Symposium on Solid Modeling and Applications, 1997, pp. 155–167.
  * [18] S. Nelaturi, V. Shapiro, Configuration products and quotients in geometric modeling, Computer-Aided Design 43 (7) (2011) 781–794.
  * [19] M. Lysenko, S. Nelaturi, V. Shapiro, Group morphology with convolution algebras, in: Proceedings of the 14th ACM Symposium on Solid and Physical Modeling, 2010, pp. 11–22.
  * [20] M. Behandish, S. Nelaturi, J. de Kleer, Automated process planning for hybrid manufacturing, Computer-Aided Design 102 (2018) 115–127.
  * [21] M. Behandish, A. M. Mirzendehdel, S. Nelaturi, A classification of topological discrepancies in additive manufacturing, Computer-Aided Design 115 (2019) 206–217.
  * [22] S. Nelaturi, M. Behandish, A. M. Mirzendehdel, J. de Kleer, Automatic support removal for additive manufacturing post processing, Computer-Aided Design 115 (2019) 135–146.
  * [23] M. Behandish, Analytic methods for geometric modeling.
  * [24] F. Stöckli, K. Shea, Topology optimization of rigid-body systems considering collision avoidance, Journal of Mechanical Design 142 (8) (2020) 081705\.
  * [25] N. S. Iyer, A. M. Mirzendehdel, S. Raghavan, Y. Jiao, E. Ulu, M. Behandish, S. Nelaturi, D. M. Robinson, PATO: Producibility-aware topology optimization using deep learning for metal additive manufacturing, arXiv preprint arXiv:2112.04552.
  * [26] A. A. G. Requicha, Representations for rigid solids: Theory, methods, and systems, ACM Computing Surveys 12 (4) (1980) 437–464.
  * [27] T. Lozano-Perez, Spatial planning: A configuration space approach, in: Autonomous Robot Vehicles, Springer, 1990, pp. 259–271.
  * [28] H. T. Ilieş, V. Shapiro, Equivalence classes for shape synthesis of moving mechanical parts, J. Comput. Inf. Sci. Eng. 4 (1) (2004) 20–27.
  * [29] K. Svanberg, The method of moving asymptotes—a new method for structural optimization, International Journal for Numerical Methods in Engineering 24 (2) (1987) 359–373.
  * [30] C. B. Morris, A. M. Mirzendehdel, M. Behandish, Topology optimization with locally evaluable complement space connectivity, in: Proceedings of the 2021 International Design Engineering Technical Conferences and Computers and Information in Engineering Conference (IDETC/CIE’2021), Vol. 85383, American Society of Mechanical Engineers (ASME), 2021, p. V03AT03A030.

##  Appendix A Shape Discretization Schemes

To enable finite representation of solids and digital computation of collision
measures, let us approximate the indicator function using a finite linear
combination:

$\rho_{S}({\mathbf{x}})=\sum_{i=1}^{n}\rho_{S,i}\phi_{i}({\mathbf{x}}),\quad\rho_{S,i}\in[0,\infty),~{}\phi_{i}:{\mathds{R}}^{\mathrm{d}}\to{\mathds{R}},$
(48)

where $\phi_{i}({\mathbf{x}})$ are local basis functions associated with a
spatial tessellation and the coefficients $\rho_{S,i}\geq 0$ are the design
variables representing an equivalence class of solids with the same
$\mathrm{d}-$measure properties, for $i=1,2,\ldots,n$. To approximately
compute the integrals in (21) and (22), we consider an asymmetric
discretization strategy in which we discretize the stationary solid (i.e., the
one to which the frame of reference is attached) via a finite volume (FV)
scheme, while we discretize the moving solid (i.e., the one whose motion is
observed) via a finite sample (FS) scheme.

For both schemes, let us immerse the solid in a uniform Cartesian grid of edge
length $\epsilon>0$, which is small enough to resolve the minimum features of
the solid. Let ${\mathbf{x}}_{i}\in{\mathds{R}}^{\mathrm{d}}$ be the
coordinates of the i${}^{\text{th}}$ vertex (i.e., $0-$cell) on the grid. Let
$C_{i}(\epsilon)\subset{\mathds{R}}^{\mathrm{d}}$ stand for the dual
$\mathrm{d}-$cell (e.g., square pixels in 2D and cubic voxels in 3D):

$C_{i}(\epsilon)\triangleq\big{[}{\mathbf{x}}_{i}-\nicefrac{{\epsilon}}{{2}},{\mathbf{x}}_{i}+\nicefrac{{\epsilon}}{{2}}\big{]}^{\mathrm{d}}$
(49)

To capture the $\mathrm{d}-$integral quantities accurately, let us associate
to the i${}^{\text{th}}$ grid vertex/cell the $\mathrm{d}-$measure of the
solid confined to $C_{i}(\epsilon)$:

$\rho_{S,i}\triangleq\mu^{\mathrm{d}}[S\cap
C_{i}(\epsilon)]=\int_{C_{i}(\epsilon)}{\mathbf{1}}_{S}({\mathbf{x}})~{}d\mu^{\mathrm{d}}[{\mathbf{x}}].$
(50)

For the FV scheme, the basis functions are selected as unit
$\mathrm{d}-$measure pulse functions, i.e., indicator functions of the cells
$C_{i}(\epsilon)$ normalized by the cell $\mathrm{d}-$measure
$\mu^{\mathrm{d}}[C_{i}(\epsilon)]=\epsilon^{\mathrm{d}}$:

$\phi^{\mathsf{FV}}_{i}({\mathbf{x}})\triangleq\frac{{\mathbf{1}}_{C_{i}(\epsilon)}({\mathbf{x}})}{\mu^{\mathrm{d}}[C_{i}(\epsilon)]}=\frac{1}{\epsilon^{\mathrm{d}}}{\mathbf{1}}_{C(\epsilon)}({\mathbf{x}}-{\mathbf{x}}_{i}),$
(51)

where
$C(\epsilon)\triangleq[+\nicefrac{{\epsilon}}{{2}},-\nicefrac{{\epsilon}}{{2}}]$.
Using these basis functions in the discrete formula in (48) gives a piecewise-
constant approximation of the indicator function over the dual grid.

For the FS scheme, on the other hand, the basis functions are selected as unit
impulse functions, i.e., shifted Dirac delta functions that “lump” a unit
$\mathrm{d}-$measure at each cell center:

$\phi^{\mathsf{FS}}_{i}({\mathbf{x}})\triangleq{\bm{\updelta}}^{\mathrm{d}}({\mathbf{x}}-{\mathbf{x}}_{i})\cong\lim_{\epsilon\to
0^{+}}\frac{1}{\epsilon^{\mathrm{d}}}{\mathbf{1}}_{C(\epsilon)}({\mathbf{x}}-{\mathbf{x}}_{i}).$
(52)

Note that $\phi^{\mathsf{FS}}_{i}({\mathbf{x}})=\lim_{\epsilon\to
0^{+}}\phi^{\mathsf{FV}}_{i}({\mathbf{x}})$. Both basis functions have a unit
$\mathrm{d}-$measure, by construction. Using these basis functions in the
discrete formula in (48) gives a lumped-measure approximation of the indicator
function over the primal grid.

Note that both FV and FS schemes are constructed such that the
$\mathrm{d}-$measures of the solid over every subset of the $\mathrm{d}-$cells
(e.g., digitized areas in 2D and volumes in 3D) are captured exactly for the
finite $\epsilon>0$, while other inexact integral properties converge to their
exact values as $\epsilon\to 0^{+}$.

Using an FV scheme to approximate the indicator functions of stationary solids
$S_{1}\subseteq\Omega_{1}$ and $S_{2}\subseteq\Omega_{2}$, we obtain:

$\displaystyle\rho^{\mathsf{FV}}_{S_{1}}({\mathbf{x}})=\sum_{i_{1}=1}^{n_{1}}\rho_{S_{1},i_{1}}\phi^{\mathsf{FV}}_{i}({\mathbf{x}}))=\frac{1}{\epsilon^{\mathrm{d}}}\sum_{i_{1}=1}^{n_{1}}\rho_{S_{1},i_{1}}{\mathbf{1}}_{C(\epsilon)}({\mathbf{x}}-{\mathbf{x}}_{i_{1}}),$
(53)
$\displaystyle\rho^{\mathsf{FV}}_{S_{2}}({\mathbf{x}})=\sum_{i_{2}=1}^{n_{2}}\rho_{S_{2},i_{2}}\phi^{\mathsf{FV}}_{i}({\mathbf{x}})=\frac{1}{\epsilon^{\mathrm{d}}}\sum_{i_{1}=1}^{n_{1}}\rho_{S_{2},i_{2}}{\mathbf{1}}_{C(\epsilon)}({\mathbf{x}}-{\mathbf{x}}_{i_{2}}),$
(54)

where $\rho_{S_{1},i_{1}},\rho_{S_{2},i_{2}}\in[0,1]$ are
$\mathrm{d}-$measures of intersections $S_{1}\cap C_{i_{1}}(\epsilon)$ and
$S_{2}\cap C_{i_{2}}(\epsilon)$, respectively., i.e., the portions of each
solid trapped inside a given cell.

Using an FS scheme to approximate the indicator functions of relatively moving
solids $S_{1,2}(t)=\tau_{2,1}(t)S_{1}$ and $S_{2,1}(t)=\tau_{1,2}(t)S_{2}$
defined in (10) and (11), we obtain:

$\displaystyle\rho^{\mathsf{FS}}_{S_{1,2}(t)}({\mathbf{x}})$
$\displaystyle=\rho^{\mathsf{FS}}_{S_{1}}(\tau_{1,2}(t){\mathbf{x}})$ (55)
$\displaystyle=\sum_{i_{1}=1}^{n_{1}}\rho_{S_{1},i_{1}}{\bm{\updelta}}^{\mathrm{d}}\left(\tau_{1,2}(t){\mathbf{x}}-{\mathbf{x}}_{i_{1}}\right),$
(56) $\displaystyle\rho^{\mathsf{FS}}_{S_{2,1}(t)}({\mathbf{x}})$
$\displaystyle=\rho^{\mathsf{FS}}_{S_{2}}(\tau_{2,1}(t){\mathbf{x}})$ (57)
$\displaystyle=\sum_{i_{2}=1}^{n_{2}}\rho_{S_{2},i_{2}}{\bm{\updelta}}^{\mathrm{d}}\left(\tau_{2,1}(t){\mathbf{x}}-{\mathbf{x}}_{i_{2}}\right),$
(58)

Substituting the indicator functions of the stationary and moving solids in
(21) and (22) with their $\mathrm{d}-$measure-preserving approximations in
(53) through (58), respectively, and rearranging sums and integrals, we
obtain:

$\displaystyle\mathcal{g}_{2,1}$ $\displaystyle=\lim_{\epsilon\to
0^{+}}\int_{0}^{1}\int_{\Omega_{1}}\rho^{\mathsf{FS}}_{S_{2,1}(t)}({\mathbf{x}})\rho^{\mathsf{FV}}_{S_{1}}({\mathbf{x}})~{}d\mu^{\mathrm{d}}[{\mathbf{x}}]~{}dt$
$\displaystyle=\lim_{\epsilon\to
0^{+}}\sum_{i_{1}=1}^{n_{1}}\sum_{i_{2}=1}^{n_{2}}\rho_{S_{1},i_{1}}\rho_{S_{2},i_{2}}w^{1,2}_{i_{1},i_{2}},$
(59) $\displaystyle\mathcal{g}_{1,2}$ $\displaystyle=\lim_{\epsilon\to
0^{+}}\int_{0}^{1}\int_{\Omega_{2}}\rho^{\mathsf{FS}}_{S_{1,2}(t)}({\mathbf{x}})\rho^{\mathsf{FV}}_{S_{2}}({\mathbf{x}})~{}d\mu^{\mathrm{d}}[{\mathbf{x}}]~{}dt$
$\displaystyle=\lim_{\epsilon\to
0^{+}}\sum_{i_{1}=1}^{n_{1}}\sum_{i_{2}=1}^{n_{2}}\rho_{S_{1},i_{1}}\rho_{S_{2},i_{2}}w^{2,1}_{i_{2},i_{1}},$
(60)

where the weights $w^{1,2}_{i_{1},i_{2}}$ and $w^{2,1}_{i_{2},i_{1}}$ are
defined by:

$\displaystyle w^{1,2}_{i_{1},i_{2}}$
$\displaystyle\triangleq\frac{1}{\epsilon^{\mathrm{d}}}\int_{0}^{1}\int_{\Omega_{1}}{\bm{\updelta}}^{\mathrm{d}}\left(\tau_{2,1}(t){\mathbf{x}}-{\mathbf{x}}_{i_{2}}\right){\mathbf{1}}_{C_{i_{1}}}({\mathbf{x}})~{}d\mu^{\mathrm{d}}[{\mathbf{x}}]~{}dt,$
$\displaystyle w^{2,1}_{i_{2},i_{1}}$
$\displaystyle\triangleq\frac{1}{\epsilon^{\mathrm{d}}}\int_{0}^{1}\int_{\Omega_{2}}{\bm{\updelta}}^{\mathrm{d}}\left(\tau_{1,2}(t){\mathbf{x}}-{\mathbf{x}}_{i_{1}}\right){\mathbf{1}}_{C_{i_{2}}}({\mathbf{x}})~{}d\mu^{\mathrm{d}}[{\mathbf{x}}]~{}dt.$

The $\delta^{\mathrm{d}}-$function has a “sifting” property; meaning that it
turns $\mathrm{d}-$integrals into a reading of the integrand at the impulse
centers (where the $\delta^{\mathrm{d}}-$function’s input is zero):

$\displaystyle w^{1,2}_{i_{1},i_{2}}$
$\displaystyle=\frac{1}{\epsilon^{\mathrm{d}}}\int_{0}^{1}{\mathbf{1}}_{C_{i_{1}}}\left(\tau_{1,2}(t){\mathbf{x}}_{i_{2}}\right)~{}dt,$
(61) $\displaystyle w^{2,1}_{i_{2},i_{1}}$
$\displaystyle=\frac{1}{\epsilon^{\mathrm{d}}}\int_{0}^{1}{\mathbf{1}}_{C_{i_{2}}}\left(\tau_{2,1}(t){\mathbf{x}}_{i_{1}}\right)~{}dt.$
(62)

To pre-compute the grid correlations (i.e., weights) in (61) and (62), the
time integral can be discretized using a simple Riemann sum approximation:

$\displaystyle w^{1,2}_{i_{1},i_{2}}$ $\displaystyle=\lim_{\delta\to
0^{+}}\frac{\delta}{\epsilon^{\mathrm{d}}}\sum_{k=1}^{K}{\mathbf{1}}_{C_{i_{1}}}\left(\tau_{1,2}(t_{k}){\mathbf{x}}_{i_{2}}\right),$
(63) $\displaystyle w^{2,1}_{i_{2},i_{1}}$ $\displaystyle=\lim_{\delta\to
0^{+}}\frac{\delta}{\epsilon^{\mathrm{d}}}\sum_{k=1}^{K}{\mathbf{1}}_{C_{i_{2}}}\left(\tau_{2,1}(t_{k}){\mathbf{x}}_{i_{1}}\right),$
(64)

where $t_{k}\triangleq(k+\nicefrac{{1}}{{2}})\delta$ for $k=0,1,\ldots,K-1$ to
uniformly discretize the time period $[0,1]$ with time-steps of
$\delta=\nicefrac{{1}}{{K}}$, small enough to capture the minimum features of
the motion trajectories. Computing each sum takes $O(K)$ operations to
evaluate and apply the motions $\tau_{1,2}(t_{k})\in M_{1,2}$ and
$\tau_{2,1}(t_{k})\in M_{2,1}$ to each of the grid vertices
${\mathbf{x}}_{i_{1}}\in\Omega_{1}$ and ${\mathbf{x}}_{i_{2}}\in\Omega_{2}$,
and testing the membership of the displaced points against $C_{i_{1}}$ and
$C_{i_{2}}$, respectively. For all time-steps and grid points, this takes
$O(n_{1}n_{2}K)$ operations, which is carried out once in a pre-processing
step and can be perfectly parallelized on the CPU/GPU.

##  Appendix B Collision Sensitivity Analysis

The additive/subtractive topological sensitivity fields (TSFs) (denoted by
superscripts $+$/$-$, respectively) for the global collision measures
$\mathcal{g}_{2,1}$ and $\mathcal{g}_{1,2}$ in (21) and (22), respectively,
are defined in terms of the changes in measures with respect to infinitesimal
inclusion/exclusion in the shapes, for a given relative motion.

  * •

With respect to changes in the stationary shapes, the TSFs are defined by:

$\displaystyle{\mathcal{D}}^{+}_{1}\mathcal{g}_{2,1}({\mathbf{x}})$
$\displaystyle\triangleq\lim_{\epsilon\to
0^{+}}\frac{\mathcal{g}_{2,1^{+}}({\mathbf{x}})-\mathcal{g}_{2,1}}{\epsilon^{\mathrm{d}}},$
(65) $\displaystyle{\mathcal{D}}^{-}_{1}\mathcal{g}_{2,1}({\mathbf{x}})$
$\displaystyle\triangleq\lim_{\epsilon\to
0^{+}}\frac{\mathcal{g}_{2,1}-\mathcal{g}_{2,1^{-}}({\mathbf{x}})}{\epsilon^{\mathrm{d}}},$
(66)

where $\mathcal{g}_{2,1^{+}}({\mathbf{x}})$ and
$\mathcal{g}_{2,1^{-}}({\mathbf{x}})$ are the global collision measures for
the pairs
$\big{(}S_{1}\cup[{\mathbf{x}}-\nicefrac{{\epsilon}}{{2}},{\mathbf{x}}+\nicefrac{{\epsilon}}{{2}}]^{\mathrm{d}},S_{2}\big{)}$
and
$\big{(}S_{1}-({\mathbf{x}}-\nicefrac{{\epsilon}}{{2}},{\mathbf{x}}+\nicefrac{{\epsilon}}{{2}})^{\mathrm{d}},S_{2}\big{)}$,
respectively. Similarly:

$\displaystyle{\mathcal{D}}^{+}_{2}\mathcal{g}_{1,2}({\mathbf{x}})$
$\displaystyle\triangleq\lim_{\epsilon\to
0^{+}}\frac{\mathcal{g}_{1,2^{+}}({\mathbf{x}})-\mathcal{g}_{1,2}}{\epsilon^{\mathrm{d}}},$
$\displaystyle{\mathcal{D}}^{-}_{2}\mathcal{g}_{1,2}({\mathbf{x}})$
$\displaystyle\triangleq\lim_{\epsilon\to
0^{+}}\frac{\mathcal{g}_{1,2}-\mathcal{g}_{1,2^{-}}({\mathbf{x}})}{\epsilon^{\mathrm{d}}},$

where $\mathcal{g}_{1,2^{+}}({\mathbf{x}})$ and
$\mathcal{g}_{1,2^{-}}({\mathbf{x}})$ are the global collision measures for
the pairs
$\big{(}S_{2}\cup[{\mathbf{x}}-\nicefrac{{\epsilon}}{{2}},{\mathbf{x}}+\nicefrac{{\epsilon}}{{2}}]^{\mathrm{d}},S_{1}\big{)}$
and
$\big{(}S_{2}-({\mathbf{x}}-\nicefrac{{\epsilon}}{{2}},{\mathbf{x}}+\nicefrac{{\epsilon}}{{2}})^{\mathrm{d}},S_{1}\big{)}$,
respectively.

  * •

With respect to changes in the moving shapes, the TSFs are defined by:

$\displaystyle{\mathcal{D}}^{+}_{2}\mathcal{g}_{2,1}({\mathbf{x}})$
$\displaystyle\triangleq\lim_{\epsilon\to
0^{+}}\frac{\mathcal{g}_{2^{+},1}({\mathbf{x}})-\mathcal{g}_{2,1}}{\epsilon^{\mathrm{d}}},$
$\displaystyle{\mathcal{D}}^{-}_{2}\mathcal{g}_{2,1}({\mathbf{x}})$
$\displaystyle\triangleq\lim_{\epsilon\to
0^{+}}\frac{\mathcal{g}_{2,1}-\mathcal{g}_{2^{-},1}({\mathbf{x}})}{\epsilon^{\mathrm{d}}},$

where $\mathcal{g}_{2^{+},1}({\mathbf{x}})$ and
$\mathcal{g}_{2^{-},1}({\mathbf{x}})$ are the global collision measures for
the pairs
$\big{(}S_{1},S_{2}\cup[{\mathbf{x}}-\nicefrac{{\epsilon}}{{2}},{\mathbf{x}}+\nicefrac{{\epsilon}}{{2}}]^{\mathrm{d}}\big{)}$
and
$\big{(}S_{1},S_{2}-({\mathbf{x}}-\nicefrac{{\epsilon}}{{2}},{\mathbf{x}}+\nicefrac{{\epsilon}}{{2}})^{\mathrm{d}}\big{)}$,
respectively. Similarly:

$\displaystyle{\mathcal{D}}^{+}_{1}\mathcal{g}_{1,2}({\mathbf{x}})$
$\displaystyle\triangleq\lim_{\epsilon\to
0^{+}}\frac{\mathcal{g}_{1^{+},2}({\mathbf{x}})-\mathcal{g}_{1,2}}{\epsilon^{\mathrm{d}}},$
$\displaystyle{\mathcal{D}}^{-}_{1}\mathcal{g}_{1,2}({\mathbf{x}})$
$\displaystyle\triangleq\lim_{\epsilon\to
0^{+}}\frac{\mathcal{g}_{1,2}-\mathcal{g}_{1^{-},2}({\mathbf{x}})}{\epsilon^{\mathrm{d}}},$

where $\mathcal{g}_{1^{+},2}({\mathbf{x}})$ and
$\mathcal{g}_{1^{-},2}({\mathbf{x}})$ are the global collision measures for
the pairs
$\big{(}S_{2},S_{1}\cup[{\mathbf{x}}-\nicefrac{{\epsilon}}{{2}},{\mathbf{x}}+\nicefrac{{\epsilon}}{{2}}]^{\mathrm{d}}\big{)}$
and
$\big{(}S_{2},S_{1}-({\mathbf{x}}-\nicefrac{{\epsilon}}{{2}},{\mathbf{x}}+\nicefrac{{\epsilon}}{{2}})^{\mathrm{d}}\big{)}$,
respectively.

Hereafter, we develop the relationships for the first two TSFs
${\mathcal{D}}^{+}_{1}\mathcal{g}_{2,1}({\mathbf{x}})$ and
${\mathcal{D}}^{-}_{1}\mathcal{g}_{2,1}({\mathbf{x}})$ defined in (65) and
(66), respectively, for sensitivity analysis with respect to an infinitesimal
inclusion/exclusion in the stationary solid $S_{1}$. The same relationships
hold for the remaining six TSFs due to the symmetry in formulation with
respect to index assignment and choice of frame of reference.

It is easy to verify that the additive/subtractive TSFs vanish when the query
point is inside/outside, respectively, of the solid with respect to which the
TSFs are computed, because the solid does not change due to union/inclusion
with a small neighborhood:

$\displaystyle\exists\epsilon>0~{}\text{s.t.}~{}\mathcal{g}_{2,1^{+}}({\mathbf{x}})=\mathcal{g}_{2,1}$
$\displaystyle~{}\Rightarrow~{}{\mathcal{D}}^{+}_{1}\mathcal{g}_{2,1}({\mathbf{x}})=0,\quad\text{if}~{}{\mathbf{x}}\in\mathsf{i}S_{1},$
$\displaystyle\exists\epsilon>0~{}\text{s.t.}~{}\mathcal{g}_{2,1^{-}}({\mathbf{x}})=\mathcal{g}_{2,1}$
$\displaystyle~{}\Rightarrow~{}{\mathcal{D}}^{-}_{1}\mathcal{g}_{2,1}({\mathbf{x}})=0,\quad\text{if}~{}{\mathbf{x}}\in\mathsf{e}S_{1},$

noting that a solid $S_{1}$ is a closed regular set, i.e., contains both its
interior (an open set) $\mathsf{i}S_{1}=S_{1}-\partial S_{1}$ and boundary
$\partial S_{1}$, while the exterior is the complement (also an open set):
$\mathsf{e}S_{1}=\overline{S}_{1}={\mathds{R}}^{\mathrm{d}}-S_{1}$. The
interior, boundary, and exterior partition the $\mathrm{d}-$space (in/on/out
classification).

The collision integral in (21) is an additive property, i.e., the change in
the integral due to the change in shape can be computed by applying the
integral to the change itself:

$\displaystyle{\mathcal{D}}^{+}_{1}\mathcal{g}_{2,1}({\mathbf{x}})$
$\displaystyle=\lim_{\epsilon\to
0^{+}}\frac{1}{\epsilon^{\mathrm{d}}}\int_{0}^{1}\int_{C({\mathbf{x}},\epsilon)-S_{1}}{\mathbf{1}}_{S_{2,1}(t)}({\mathbf{x}}^{\prime})~{}d\mu^{\mathrm{d}}[{\mathbf{x}}^{\prime}]~{}dt,$
(67) $\displaystyle{\mathcal{D}}^{-}_{1}\mathcal{g}_{2,1}({\mathbf{x}})$
$\displaystyle=\lim_{\epsilon\to
0^{+}}\frac{1}{\epsilon^{\mathrm{d}}}\int_{0}^{1}\int_{C({\mathbf{x}},\epsilon)\cap
S_{1}}{\mathbf{1}}_{S_{2,1}(t)}({\mathbf{x}}^{\prime})~{}d\mu^{\mathrm{d}}[{\mathbf{x}}^{\prime}]~{}dt,$
(68)

where
$C({\mathbf{x}},\epsilon)\triangleq({\mathbf{x}}-\nicefrac{{\epsilon}}{{2}},{\mathbf{x}}+\nicefrac{{\epsilon}}{{2}})^{\mathrm{d}}$.

For interior points ${\mathbf{x}}_{\mathsf{i}}\in\mathsf{i}S_{1}$, we have
$C({\mathbf{x}}_{\mathsf{i}},\epsilon)\subset S_{1}$ for small enough
$\epsilon>0$, hence $C({\mathbf{x}}_{\mathsf{i}},\epsilon)-S_{1}=\emptyset$
and the integral in (67) vanishes, as expected. But
$C({\mathbf{x}}_{\mathsf{i}},\epsilon)\cap
S_{1}=C({\mathbf{x}}_{\mathsf{i}},\epsilon)$ and the integral in (68) reduces
to:

${\mathcal{D}}^{-}_{1}\mathcal{g}_{2,1}({\mathbf{x}}_{\mathsf{i}})=\lim_{\epsilon\to
0^{+}}\int_{0}^{1}\int_{C({\mathbf{x}}_{\mathsf{i}},\epsilon)}\frac{1}{\epsilon^{\mathrm{d}}}{\mathbf{1}}_{S_{2,1}(t)}({\mathbf{x}}^{\prime})~{}d\mu^{\mathrm{d}}[{\mathbf{x}}^{\prime}]~{}dt.$
(69)

If, additionally, $\tau_{2,1}(t){\mathbf{x}}_{\mathsf{i}}\in\mathsf{i}S_{2}$,
i.e., ${\mathbf{x}}_{\mathsf{i}}\in\tau_{1,2}(t)\mathsf{i}S_{2}$, then
$C({\mathbf{x}}_{\mathsf{i}},\epsilon)\subset\tau_{1,2}(t)\mathsf{i}S_{2}=\mathsf{i}S_{2,1}(t)$
for small enough $\epsilon>0$, hence the inner integral equals
$\epsilon^{\mathrm{d}}/\epsilon^{\mathrm{d}}=1$.

For exterior points ${\mathbf{x}}_{\mathsf{e}}\in\mathsf{e}S_{1}$, we have
$C({\mathbf{x}}_{\mathsf{e}},\epsilon)\subset S_{1}$ for small enough
$\epsilon>0$, hence $C({\mathbf{x}}_{\mathsf{e}},\epsilon)\cap
S_{1}=\emptyset$ and the integral in (68) vanishes, as expected. But
$C({\mathbf{x}}_{\mathsf{e}},\epsilon)-S_{1}=C({\mathbf{x}}_{\mathsf{e}},\epsilon)$
and the integral in (67) reduces to:

${\mathcal{D}}^{+}_{1}\mathcal{g}_{2,1}({\mathbf{x}}_{\mathsf{i}})=\lim_{\epsilon\to
0^{+}}\int_{0}^{1}\int_{C({\mathbf{x}}_{\mathsf{e}},\epsilon)}\frac{1}{\epsilon^{\mathrm{d}}}{\mathbf{1}}_{S_{2,1}(t)}({\mathbf{x}}^{\prime})~{}d\mu^{\mathrm{d}}[{\mathbf{x}}^{\prime}]~{}dt.$
(70)

If, additionally, $\tau_{2,1}(t){\mathbf{x}}_{\mathsf{e}}\in\mathsf{i}S_{2}$,
i.e., ${\mathbf{x}}_{\mathsf{e}}\in\tau_{1,2}(t)\mathsf{i}S_{2}$, then
$C({\mathbf{x}}_{\mathsf{e}},\epsilon)\subset\tau_{1,2}(t)\mathsf{i}S_{2}=\mathsf{i}S_{2,1}(t)$
for small enough $\epsilon>0$, hence the inner integral equals
$\epsilon^{\mathrm{d}}/\epsilon^{\mathrm{d}}=1$.

To summarize, if either ${\mathbf{x}}\in\mathsf{i}S_{1}$ or
${\mathbf{x}}\in\mathsf{e}S_{1}$, then:

  * •

If $\tau_{2,1}(t){\mathbf{x}}\in\mathsf{i}S_{2}$, then the inner integral is
1.

  * •

If $\tau_{2,1}(t){\mathbf{x}}\in\mathsf{e}S_{2}$, then the inner integral is
0.

The case for query points that fall on the boundaries of either shape is
slightly more complicated.

Assuming the curve $\tau_{2,1}(t){\mathbf{x}}$ intersects the boundary
$\partial S_{2}$ over singular points (as opposed to continuous curve
segments) the outer integral measures the duration of time over which the
trajectory collides with the solid $S_{2}$ as $\epsilon\to 0^{+}$. Hence, the
limit converges and both additive and subtractive TSFs exists.444If the curve
remains tangent to the boundary over a continuous curve segment, the duration
of time over which this happens is weighted depending on the local
neighborhood shape in the limit $\epsilon\to 0^{+}$, but the fact that TSF is
convergent does not change.

Note also that the the formulas for additive and subtractive TSFs are
identical, hence they converge to the same formula as
$\|{\mathbf{x}}_{\mathsf{e}}-{\mathbf{x}}_{\mathsf{i}}\|_{2}\to 0^{+}$,
meaning that the interior and exterior query points converge to a common
meeting point ${\mathbf{x}}\in\partial S_{1}$ on the boundary. However, the
value of ${\mathcal{D}}^{+}_{1}\mathcal{g}_{2,1}({\mathbf{x}})$ and
${\mathcal{D}}^{-}_{1}\mathcal{g}_{2,1}({\mathbf{x}})$ when the query point is
precisely on the boundary ${\mathbf{x}}\in\partial S_{1}$ may be different,
depending on the shape of the complementary partial neighborhoods
$C({\mathbf{x}},\epsilon)-S_{1}$ and $C({\mathbf{x}},\epsilon)\cap S_{1}$. In
other words, even though the TSFs exist and have limit relationships, they are
discontinuous as the query point passes through the boundary.

Next, let us consider the discretization schemes presented in 3.2 and A and
how the TSFs relate to the discrete derivatives given in (31) and (32) as
$\epsilon\to 0^{+}$. Remember that the global collision measures can be
approximated by (27) and (28) in the following sense:

$\displaystyle\mathcal{g}_{2,1}$ $\displaystyle=\lim_{\epsilon\to
0^{+}}\sum_{i_{1}=1}^{n_{1}}\sum_{i_{2}=1}^{n_{2}}\rho_{S_{1},i_{1}}\rho_{S_{2},i_{2}}w^{1,2}_{i_{1},i_{2}},$
(71) $\displaystyle\mathcal{g}_{1,2}$ $\displaystyle=\lim_{\epsilon\to
0^{+}}\sum_{i_{1}=1}^{n_{1}}\sum_{i_{2}=1}^{n_{2}}\rho_{S_{1},i_{1}}\rho_{S_{2},i_{2}}w^{2,1}_{i_{2},i_{1}},$
(72)

where $\rho_{S_{1},i_{1}},\rho_{S_{2},i_{2}}\geq 0$ were defined by (50).
Applying the additive/subtractive TSF operators to both sides of (71), which
we explored above, yields:

${\mathcal{D}}^{\pm}_{1}\mathcal{g}_{2,1}({\mathbf{x}})=\lim_{\epsilon\to
0^{+}}\sum_{i_{1}=1}^{n_{1}}\sum_{i_{2}=1}^{n_{2}}\left({\mathcal{D}}^{\pm}_{1}\rho_{S_{1},i_{1}}({\mathbf{x}})\right)\rho_{S_{2},i_{2}}w^{1,2}_{i_{1},i_{2}},$
(73)

noting that the infinitesimal inclusion/exclusion are applied to $S_{1}$ in
this case, thus they affect $\rho_{S_{1},i_{1}}$ but not $\rho_{S_{2},i_{2}}$.

Substituting for (50) and, once again, exploiting the additivity of the
integral definition, we obtain:

$\displaystyle{\mathcal{D}}^{+}_{1}\rho_{S_{1},i_{1}}({\mathbf{x}})$
$\displaystyle=\lim_{\epsilon^{\prime}\to
0^{+}}\frac{1}{(\epsilon^{\prime})^{\mathrm{d}}}\int_{C({\mathbf{x}}_{i_{1}},\epsilon)-S_{1}}{\mathbf{1}}_{C({\mathbf{x}},\epsilon^{\prime})}({\mathbf{x}}^{\prime})~{}d\mu^{\mathrm{d}}[{\mathbf{x}}^{\prime}],$
$\displaystyle{\mathcal{D}}^{-}_{1}\rho_{S_{1},i_{1}}({\mathbf{x}})$
$\displaystyle=\lim_{\epsilon^{\prime}\to
0^{+}}\frac{1}{(\epsilon^{\prime})^{\mathrm{d}}}\int_{C({\mathbf{x}}_{i_{1}},\epsilon)\cap
S_{1}}{\mathbf{1}}_{C({\mathbf{x}},\epsilon^{\prime})}({\mathbf{x}}^{\prime})~{}d\mu^{\mathrm{d}}[{\mathbf{x}}^{\prime}].$

Once again, the following two cases are simple:

  * •

If the finite grid cell is completely inside $S_{1}$, then
$C({\mathbf{x}}_{i_{1}},\epsilon)\cap S_{1}=C({\mathbf{x}}_{i_{1}},\epsilon)$,
hence:

${\mathcal{D}}^{-}_{1}\rho_{S_{1},i_{1}}({\mathbf{x}})=\lim_{\epsilon^{\prime}\to
0^{+}}\int_{C({\mathbf{x}}_{i_{1}},\epsilon)}\frac{1}{(\epsilon^{\prime})^{\mathrm{d}}}{\mathbf{1}}_{C({\mathbf{x}},\epsilon^{\prime})}({\mathbf{x}}^{\prime})~{}d\mu^{\mathrm{d}}[{\mathbf{x}}^{\prime}].$
(74)

and ${\mathcal{D}}^{+}_{1}\rho_{S_{1},i_{1}}({\mathbf{x}})=0$ because
$C({\mathbf{x}}_{i_{1}},\epsilon)-S_{1}=\emptyset$.

  * •

If the finite grid cell is completely outside $S_{1}$, then
$C({\mathbf{x}}_{i_{1}},\epsilon)-S_{1}=C({\mathbf{x}}_{i_{1}},\epsilon)$,
hence:

${\mathcal{D}}^{+}_{1}\rho_{S_{1},i_{1}}({\mathbf{x}})=\lim_{\epsilon^{\prime}\to
0^{+}}\int_{C({\mathbf{x}}_{i_{1}},\epsilon)}\frac{1}{(\epsilon^{\prime})^{\mathrm{d}}}{\mathbf{1}}_{C({\mathbf{x}},\epsilon^{\prime})}({\mathbf{x}}^{\prime})~{}d\mu^{\mathrm{d}}[{\mathbf{x}}^{\prime}].$
(75)

and ${\mathcal{D}}^{-}_{1}\rho_{S_{1},i_{1}}({\mathbf{x}})=0$ because
$C({\mathbf{x}}_{i_{1}},\epsilon)\cap S_{1}=\emptyset$.

Once again, the situation is more complicated when the grid cell cuts through
the boundary, i.e., $C({\mathbf{x}}_{i_{1}},\epsilon)\cap\partial
S_{1}\neq\emptyset$.

Note the subtle differences among the inclusion/exclusion center
${\mathbf{x}}\in{\mathds{R}}^{\mathrm{d}}$ (i.e., the point at which TSFs are
queried), the grid cell center
${\mathbf{x}}_{i_{1}}\in{\mathds{R}}^{\mathrm{d}}$ (used to discretize
$S_{1}$), the integration variable
${\mathbf{x}}^{\prime}\in{\mathds{R}}^{\mathrm{d}}$, as well as the
distinction between the grid cell size $\epsilon>0$ (finite in (74) and (75)
for a given grid, despite infinitesimal in (73)) and the inclusion/exclusion
size $\epsilon^{\prime}>0$ (infinitesimal).

Now, if we let ${\mathbf{x}}={\mathbf{x}}_{i_{1}}$, for small enough
$\epsilon^{\prime}>0$ we will have $\epsilon^{\prime}<\epsilon$ thus
$C({\mathbf{x}},\epsilon^{\prime})\subset C({\mathbf{x}},\epsilon)$ and the
identical integrals in (74) and (74) yield
$(\epsilon^{\prime})^{\mathrm{d}}/(\epsilon^{\prime})^{\mathrm{d}}=1$, thus
(73) yields:

${\mathcal{D}}^{\pm}_{1}\mathcal{g}_{2,1}({\mathbf{x}}_{i})=\lim_{\epsilon\to
0^{+}}\sum_{i_{1}=1}^{n_{1}}\sum_{i_{2}=1}^{n_{2}}\rho_{S_{2},i_{2}}w^{1,2}_{i_{1},i_{2}},$
(76)

The right-hand side is the limit of the discrete derivative
$\partial\mathcal{g}_{2,1}/\partial\rho_{S_{1},i_{1}}$ in (31), proving that
the discrete sensitivity converges to the topological sensitivity for the
fully-internal and fully-external grid cells. If the cell center is not on the
boundary, even if the cell cuts through the boundary, it will become fully-
internal or fully-external for small enough $\epsilon>0$ (i.e., resolution
refinement), thus the convergence holds. If the cell center is precisely on
the boundary, the derivation becomes slightly more involved, but follows a
similar logic.

Note also that if two solids are in collision-free contact, the TSFs of the
global collision measures queried at the boundary will be nonzero. The reason
is that a small perturbation to the shape at the contact interface may result
in collision. While this may create conservative material removal based on
discrete sensitivity analysis (e.g., via TO), such issues are inevitable due
to the loss of information as a result of discretization, e.g., ambiguity at
the boundaries where pixels (2D) or voxels (3D) approximate an uncountably
many different geometries with the same area (2D) or volume (3D). These issues
can only be ameliorated by resolution refinement.