Topology optimization with physics-informed neural networks: application
  to noninvasive detection of hidden geometries
Detecting hidden geometrical structures from surface measurements under
electromagnetic, acoustic, or mechanical loading is the goal of noninvasive
imaging techniques in medical and industrial applications. Solving the inverse
problem can be challenging due to the unknown topology and geometry, the
sparsity of the data, and the complexity of the physical laws. Physics-informed
neural networks (PINNs) have shown promise as a simple-yet-powerful tool for
problem inversion, but they have yet to be applied to general problems with a
priori unknown topology. Here, we introduce a topology optimization framework
based on PINNs that solves geometry detection problems without prior knowledge
of the number or types of shapes. We allow for arbitrary solution topology by
representing the geometry using a material density field that approaches binary
values thanks to a novel eikonal regularization. We validate our framework by
detecting the number, locations, and shapes of hidden voids and inclusions in
linear and nonlinear elastic bodies using measurements of outer surface
displacement from a single mechanical loading experiment. Our methodology opens
a pathway for PINNs to solve various engineering problems targeting geometry
optimization.
Saviz Mowlavi, Ken Kamrin
# Topology optimization with physics-informed neural networks:

###### Abstract

Detecting hidden geometrical structures from surface measurements under electromagnetic, acoustic, or mechanical loading is the goal of noninvasive imaging techniques in medical and industrial applications. Solving the inverse problem can be challenging due to the unknown topology and geometry, the sparsity of the data, and the complexity of the physical laws. Physics-informed neural networks (PINNs) have shown promise as a simple-yet-powerful tool for problem inversion, but they have yet to be applied to general problems with a priori unknown topology. Here, we introduce a topology optimization framework based on PINNs that solves geometry detection problems without prior knowledge of the number or types of shapes. We allow for arbitrary solution topology by representing the geometry using a material density field that approaches binary values thanks to a novel eikonal regularization. We validate our framework by detecting the number, locations, and shapes of hidden voids and inclusions in linear and nonlinear elastic bodies using measurements of outer surface displacement from a single mechanical loading experiment. Our methodology opens a pathway for PINNs to solve various engineering problems targeting geometry optimization.

Noninvasive detection of hidden geometries is desirable in countless applications including medical imaging and diagnosis [13], nondestructive evaluation of materials [28], and mine detection [49]. The goal is to infer the locations, and shapes of structures hidden inside a matrix from surface measurements of the response to an applied load such as a magnetic field [19; 22], an electric current [2], or a mechanical traction [9]. Identifying these internal boundaries from the measured data constitutes a challenging inverse problem due to the unknown topology, the large number of parameters required to describe arbitrary geometries [27], the potential sparsity of the data, and the complexity of the underlying physical laws which usually take the form of linear or nonlinear partial differential equations (PDEs) [18].

In recent years, physics-informed neural networks (PINNs) have emerged as a robust tool for problem inversion across disciplines and over a range of model complexity [16; 33; 45]. PINNs' ability to seamlessly blend measurement data or design objectives with governing PDEs in nontrivial geometries has enabled practitioners to solve easily a range of inverse problems involving identification or design of unknown properties in fields ranging from mechanics to optics and medicine [11; 25; 36; 42; 46; 50]. Encouraged by these early successes, we introduce in this paper a general topology optimization (TO) framework to solve noninvasive geometry detection problems using PINNs, leveraging both the measurements and the governing PDEs. Building on the strength of PINNs, our approach is straightforward to implement regardless of the complexity of the physical model, produces accurate results using measurement data from a single experiment, and does not require a training dataset. To the best of our knowledge, the present work marks the first time that PINNs have been applied to problems involving _a priori_ unknown topology and geometry.

Classical approaches to solve geometry identification inverse problems tend to be complex as they combine traditional numerical solvers such as the finite-element or boundary-element method, adjoint techniques to evaluate the sensitivity of the error residual with respect to the shape or the topology, and gradient descent-based optimization algorithms to update the geometry at every iteration [4; 7; 18; 20; 24]. Furthermore, these techniques call for carefully chosen regularization schemes and do not always yield satisfactory results, especially in the presence of sparse measurements acquired by only one or a few sets of experiments. For example, in the mechanical loading case where voids and inclusions in an elastic body are to be identified from measurements of surface displacement in response to a prescribed traction, past studies limit themselves to simple shapes like squares and circles or fail to find the right number of shapes [5; 9; 34; 39; 40]. Attempts to apply PINNs to geometry detection are in their infancy, so far restricted to cases where the number and shape-type of hidden voids or inclusions in an elastic structure are provided in advance [56].

Our PINN-based TO framework does not require any prior knowledge on the number and types of shapes. We allow for arbitrary solution topology by representing the geometry using a material density field equal to 0 in one phase and 1 in the other. The material density is parameterized through a neural network, which needs to be regularized in order to push the material density towards 0 or 1 values. Thus, one key ingredient in our framework is a novel eikonal regularization, inspired from fast-marching level-set methods [1; 43] and neural signed distance functions [23], that promotes a constant thickness of the interface region where the material density transitions between 0 and 1, leading to well-defined boundaries throughout the domain. This eikonal regularization en
ters as an additional term in the standard PINN loss, which is then used to train the neural networks underlying the material density and physical quantities to yield a solution to the geometry detection problem. As an illustration, we apply our framework to cases involving an elastic body under mechanical loading, and discover the topology, locations, and shapes of hidden structures for a variety of geometries and materials.

## Results

### Problem formulation

We consider noninvasive geometry detection problems of the following form. Suppose we have a continuous body \(\mathcal{B}\) containing an unknown number of hidden voids or inclusions, with unknown shapes and at unknown locations within the body. The material properties are assumed to be known and homogeneous within the body and the inclusions. We then apply a certain type of loading (e.g. mechanical, thermal, acoustic, etc) on the body's external boundary \(\partial\mathcal{B}^{\mathrm{ext}}\), which produces a response within the body that can be described by a set of \(n\) physical quantities (e.g. displacements, stresses, temperature, etc). These physical quantities can be lumped into a vector field \(\mathbf{\psi}:\mathcal{B}\rightarrow\mathbb{R}^{n}\) and satisfy a known set of governing PDEs. The goal of the inverse problem is to identify the number, locations, and shapes of the voids or inclusions based on measurements along \(\partial\mathcal{B}^{\mathrm{ext}}\) of some of the physical quantities contained in \(\mathbf{\psi}\).

As a concrete example, we consider two prototypical plane-strain elasticity inverse problems. In the first case, a square elastic matrix with hidden voids or inclusions is pulled by a uniform traction \(P_{o}\) on two sides (Fig. 1a). The goal of the inverse problem is to identify the number, locations, and shapes of the voids or inclusions using discrete measurements of the displacement of the outer boundary of the matrix. In the second case, an elastic layer on top of a hidden rigid substrate is compressed from the top by a uniform pressure \(P_{o}\), with periodic lateral boundary conditions (Fig. 1b). The goal is to identify the shape of the substrate using discrete measurements of the displacement of the top surface. For both cases, the constitutive properties of all materials are assumed to be known. We will consider two different types of constitutive laws: compressible linear elasticity, which characterizes the small deformation of any compressible elastic material, and incompressible nonlinear hyperelasticity, which models the large deformation of rubber-like materials. In the linear elastic case, there exists a unique solution to the inverse problem (see proof in Supplementary Information), making it well-suited to evaluating the accuracy of our TO framework.

Following density-based TO methods [51], we avoid any restriction on the number and shapes of hidden structures by parameterizing the geometry of the elastic body \(\mathcal{B}\) through a discrete-valued material density function \(\rho:\Omega\rightarrow\{0,1\}\), where \(\Omega\) is a global domain comprising both \(\mathcal{B}\) and the hidden voids or inclusions. The material density is defined to be equal to \(1\) in the elastic body \(\mathcal{B}\) and \(0\) in the voids or inclusions. The physical quantities \(\mathbf{\psi}\) can then be extended to the global domain \(\Omega\) by introducing an explicit \(\rho\)-dependence in their governing PDEs, leading to equations of the form

\[\mathbf{r}(\mathbf{\psi}(\mathbf{x}),\rho(\mathbf{x})) =0,\quad\mathbf{x}\in\Omega, \tag{1a}\] \[\mathbf{b}(\mathbf{\psi}(\mathbf{x})) =0,\quad\mathbf{x}\in\partial\Omega, \tag{1b}\]

with known boundary conditions defined solely on the external boundary \(\partial\Omega=\partial\mathcal{B}^{\mathrm{ext}}\). Note that the residual functions \(\mathbf{r}\) and \(\mathbf{b}\) may contain partial derivatives of \(\mathbf{\psi}\) and \(\rho\). The inverse problem is now to find the distribution of material density \(\rho\) in \(\Omega\) so that the corresponding solution for \(\mathbf{\psi}\) matches surface measurements \(\mathbf{\psi}_{i}^{m}\) at discrete locations \(\mathbf{x}_{i}\in\partial\Omega^{m}\subset\partial\Omega\), that is,

\[\mathbf{\psi}(\mathbf{x}_{i})=\mathbf{\psi}_{i}^{m},\quad\mathbf{x}_{i}\in\partial \Omega^{m}. \tag{2}\]

In practice, we might only measure select quantities in \(\mathbf{\psi}\) at some of the locations, but we do not write so explicitly to avoid overloading the notation.

For the linear elasticity problem that we consider as an example, \(\mathbf{\psi}=(\mathbf{u},\mathbf{\sigma})\) where \(\mathbf{u}(\mathbf{x})\) and \(\mathbf{\sigma}(\mathbf{x})\) are displacement and stress fields, respectively, and the governing

Figure 1: **Setup of two geometry identification problems in elastic bodies under mechanical loading.****a**, A square elastic matrix with hidden voids or inclusions is pulled by a known uniform traction on two opposite sides. The goal is to identify the number, locations, and shapes of the voids or inclusions within using measurements of the displacement occurring along the outer boundary of the matrix. **b**, An elastic layer on top of a hidden rigid substrate is compressed from the top by a uniform pressure. The goal is to identify the shape of the substrate using measurements of the displacement of the top surface.
equations comprise equilibrium relations \(\sum_{j}\partial\sigma_{ij}/\partial x_{j}=0\) and a constitutive law \(F(\mathbf{\sigma},\nabla\mathbf{u},\rho)=0\), both defined over \(\Omega\). The presence of \(\rho\) in the constitutive law specifies different material behaviors for the elastic solid phase and the void or rigid inclusion phase. The applied boundary conditions take the form \(\mathbf{u}=\bar{\mathbf{u}}\) on \(\partial\Omega_{u}\) and \(\mathbf{\sigma}\mathbf{n}=\bar{\mathbf{t}}\) on \(\partial\Omega_{t}\), where \(\partial\Omega_{u}\) and \(\partial\Omega_{t}\) are partitions of the external boundary with applied displacements and applied tractions, respectively, and \(\mathbf{n}\) is the outward unit normal. In the case of the elastic layer, the outer boundary also comprises a portion \(\partial\Omega_{p}\) with periodic boundary conditions on the displacement and traction. Finally, the requirement that the predictions for \(\mathbf{u}\) at the surface match the measurement data is expressed as \(\mathbf{u}(\mathbf{x}_{i})=\mathbf{u}_{i}^{m}\), \(\mathbf{x}_{i}\in\partial\Omega^{m}\). See Methods for a detailed formulation of the governing equations and boundary conditions for all considered cases.

Similar to density-based TO methods [51], we relax the binary constraint on the material density by allowing intermediate values of \(\rho\) between 0 and 1. This renders the problem amenable to gradient-based optimization, which underpins the PINN-based TO framework that we introduce in the next section. However, the challenge is to find an appropriate regularization mechanism that drives the optimized material distribution towards 0 and 1 rather than intermediate values devoid of physical meaning. As we will show in the discussion, common strategies employed in TO [18; 51] do not yield satisfactory results in our PINN-based framework for geometry detection problems. Thus, we have developed a novel eikonal regularization scheme inspired from level-set methods and signed distance functions, which we will describe after presenting the general framework.

### General framework

We propose a TO framework based on PINNs for solving noninvasive geometry detection problems (Fig. 2). At the core of the framework are several deep neural networks that approximate the physical quantities \(\mathbf{\psi}(\mathbf{x})\) describing the problem and the material density \(\rho(\mathbf{x})\). For the physical quantities, each neural network maps the spatial location \(\mathbf{x}=(x_{1},x_{2})\) to one of the variables in \(\mathbf{\psi}=(\psi_{1},\cdots,\psi_{n})\); this can be expressed as \(\psi_{i}=\bar{\psi}_{i}(\mathbf{x};\mathbf{\theta}_{i})\) where \(\bar{\psi}_{i}\) is the map defined by the \(i\)th neural network and its trainable parameters \(\mathbf{\theta}_{i}\) (see Methods). For the material distribution, we first define a neural network with trainable parameters \(\mathbf{\theta}_{\phi}\) that maps \(\mathbf{x}\) to a scalar variable \(\phi=\bar{\phi}(\mathbf{x};\mathbf{\theta}_{\phi})\). A sigmoid function is then applied to \(\phi\) to yield \(\rho=\text{sigmoid}(\phi/\delta)=\text{sigmoid}(\bar{\phi}(\mathbf{x};\mathbf{ \theta}_{\phi})/\delta)\), which we simply write as \(\rho=\bar{\rho}(\mathbf{x},\mathbf{\theta}_{\phi})\). This construction ensures that the material density \(\rho\) remains between 0 and 1, and \(\delta\) is a transition length scale that we will comment on later. We define the phase transition to occur at \(\rho=0.5\) so that the zero level-set of \(\phi\) delineates the boundary between the two material phases, hence \(\phi\) is hereafter referred to as a level-set function [43; 44].

We now seek the parameters \(\mathbf{\theta}_{\mathbf{\psi}}=\{\mathbf{\theta}_{1},\dots,\mathbf{\theta}_{n}\}\) and \(\mathbf{\theta}_{\phi}\) so that the neural network approximations for \(\mathbf{\psi}(\mathbf{x})\) and \(\rho(\mathbf{x})\) satisfy the governing equations (1a) and applied boundary conditions (1b) while matching the surface measurements (2). This is achieved by constructing a loss function of the form

\[\mathcal{L}(\mathbf{\theta}_{\mathbf{\psi}},\mathbf{\theta}_{\phi}) =\lambda_{\text{meas}}\mathcal{L}_{\text{meas}}(\mathbf{\theta}_{\bm {\psi}})+\lambda_{\text{gov}}\mathcal{L}_{\text{gov}}(\mathbf{\theta}_{\mathbf{\psi}},\mathbf{\theta}_{\phi})\] \[\quad+\lambda_{\text{reg}}\mathcal{L}_{\text{eik}}(\mathbf{\theta}_{ \phi}), \tag{3}\]

where \(\mathcal{L}_{\text{meas}}\) and \(\mathcal{L}_{\text{gov}}\) measure the degree to which the neural network approximations do not satisfy the measurements and governing equations, respectively, \(\mathcal{L}_{\text{eik}}\) is a crucial regularization term that drives \(\rho\) towards 0 or 1 values and that we will explain below, and the \(\lambda\)'s are scalar weights. The measurement loss takes the form

\[\mathcal{L}_{\text{meas}}(\mathbf{\theta}_{\mathbf{\psi}})=\frac{1}{|\partial\Omega^{ m}|}\sum_{\mathbf{x}_{i}\in\partial\Omega^{m}}|\bar{\mathbf{\psi}}(\mathbf{x}_{i}; \mathbf{\theta}_{\mathbf{\psi}})-\mathbf{\psi}_{i}^{m}|^{2}, \tag{4}\]

where \(|\partial\Omega^{m}|\) denotes the size of the set \(\partial\Omega^{m}\). A trivial modification of this expression is necessary in the case where only select quantities in \(\mathbf{\psi}\) are measured. The governing equations loss takes the form

\[\mathcal{L}_{\text{gov}}(\mathbf{\theta}_{\mathbf{\psi}},\mathbf{\theta}_{\phi})=\frac{1}{| \Omega^{d}|}\sum_{\mathbf{x}_{i}\in\Omega^{d}}|\mathbf{r}(\bar{\mathbf{\psi}}( \mathbf{x}_{i};\mathbf{\theta}_{\mathbf{\psi}}),\bar{\rho}(\mathbf{x}_{i};\mathbf{\theta}_ {\phi}))|^{2}, \tag{5}\]

where \(\Omega^{d}\) is a set of collocation points in \(\Omega\), and we use automatic differentiation to calculate in a mesh-free fashion the spatial derivatives contained in \(\mathbf{r}\). We design the architecture of our neural networks in such a way that they inherently satisfy the boundary conditions (see Methods, section "Detailed PINNs formulation").

Finally, the optimal parameters \(\mathbf{\theta}_{\mathbf{\psi}}^{*}\) and \(\mathbf{\theta}_{\phi}^{*}\) that solve the problem can be obtained by training the neural networks to minimize the loss (3) using stochastic gradient descent-based optimization. The corresponding physical quantities \(\bar{\mathbf{\psi}}(\mathbf{x};\mathbf{\theta}_{\mathbf{\psi}}^{*})\) will match the discrete surface measurements while satisfying the governing equations of the problem, while the corresponding material density \(\bar{\rho}(\mathbf{x};\mathbf{\theta}_{\phi}^{*})\) will reveal the number, locations, and shapes of the hidden voids or inclusions.

### Material density regularization

We now describe the key ingredient that ensures the success of our framework. As mentioned above, the main challenge is to promote the material density \(\rho(\mathbf{x})\) to converge towards 0 or 1 away from the material phase boundaries, given by the zero level-set \(\phi=0\). Moreover, we desire the thickness of the transition region along these boundaries, where \(\rho\) goes from 0 to 1, to be uniform everywhere in order to ensure consistency of physical laws across the interface (e.g. stress jumps).

To visualize what happens in the absence of regularization, consider a random instance of the neural network
\(\phi=\bar{\phi}(\mathbf{x},\mathbf{\theta}_{\phi})\) (Fig. 3a, left) and the corresponding material distribution \(\rho=\text{sigmoid}(\phi/\delta)\) with \(\delta=0.01\) (Fig. 3a, center). The sigmoid transformation ensures that \(\rho\) never drops below 0 or exceeds 1, leading to large regions corresponding to one phase or the other. However, the thickness of the transition region where \(\rho\) goes from 0 to 1 is not everywhere uniform, resulting in large zones where \(\rho\) assumes nonphysical values between 0 and 1 (Fig. 3a, center). This behavior stems from the non-uniformity of the gradient norm \(|\nabla\phi|\) along the material boundaries \(\phi=0\), with small and large values of \(|\nabla\phi|\) leading to wide and narrow transition regions, respectively (Fig. 3a, right).

We propose to regularize the material density by forcing the gradient norm \(|\nabla\phi|\) to be unity in a narrow band \(\Omega_{\text{eik}}\) of width \(w\) along the material boundaries defined by the zero level-set \(\phi=0\). In this way, \(\phi\) becomes a signed distance function to the material boundary in the narrow band, thereby constraining the gradient of \(\rho\) to be constant along the interface. To ensure that the narrow band covers the near-entirety of the transition region where \(\rho\) goes from 0 to 1, we choose \(w=10\delta\) so that \(\rho=\text{sigmoid}(\pm w/2\delta)=\text{sigmoid}(\pm 5)\simeq 0\) or 1 along the edge of the narrow band. To illustrate the effect of such regularization, we consider the previous random instance of the neural network \(\phi=\bar{\phi}(\mathbf{x},\mathbf{\theta}_{\phi})\) and enforce the constraint \(|\nabla\phi|=1\) in the narrow band \(\Omega_{\text{eik}}\) along its zero level-set (Fig. 3b, left and right). The zero level-set is kept fixed to facilitate comparison with the un-regularized case (Fig. 3a). With \(\phi\) now behaving like a signed-distance function in the narrow band, a uniform transition thickness for \(\rho\) along all material boundaries is achieved, without large regions of intermediate density values (Fig. 3b, center).

In practice, we implement this regularization into our PINN-based TO framework by including an 'eikonal' loss term \(\mathcal{L}_{\text{eik}}\) in (3), which takes the form

\[\mathcal{L}_{\text{eik}}(\mathbf{\theta}_{\phi})=\frac{1}{|\Omega_{\text{eik}}^{ d}|}\sum_{\mathbf{x}_{i}\in\Omega_{\text{eik}}^{d}}\left(|\nabla\phi(\mathbf{x}_{i} )|-1\right)^{2}, \tag{6}\]

where \(\Omega_{\text{eik}}^{d}=\{\mathbf{x}_{i}\in\Omega^{d}:|\phi(\mathbf{x}_{i})|<w /2\}\). The aim of this term is to penalize deviations away from the constraint \(|\nabla\phi|=1\) in the narrow band \(\Omega_{\text{eik}}\) of width \(w\) along the interface defined by the zero level-set \(\phi=0\). Because finding the subset of collocation points \(\mathbf{x}_{i}\) in \(\Omega^{d}\) belonging to the true narrow band of width \(w\) at every step of the training process would be too expensive, we instead relax the domain over which the constraint \(|\nabla\phi|=1\) is active by utilizing the subset \(\Omega_{\text{eik}}^{d}\) of collocation points that satisfy \(|\phi(\mathbf{x}_{i})|<w/2\). As the constraint \(|\nabla\phi|=1\) is progressively better satisfied during the training process, \(\Omega_{\text{eik}}^{d}\) will eventually overlap the true narrow band of width \(w\) along the zero level-set of \(\phi\) (Fig. 3c).

Since the constraint \(|\nabla\phi|=1\) in the narrow band takes the form of an eikonal equation, we call this approach eikonal regularization. We emphasize that in contrast to recent works training neural networks to solve the eikonal equation [23], our eikonal regularization does not force \(\phi\)

Figure 2: **TO framework for noninvasive detection of hidden geometries.** The geometry of the system, which is initially unknown, is parameterized by a material density field given through a level-set function and equal to 1 in the elastic body and 0 in the voids or inclusions. The level-set function and the physical quantities describing the problem are approximated with deep neural networks designed to inherently satisfy the applied boundary conditions. These neural networks are then trained to minimize a loss function that drives the material density and physical quantities towards satisfying the governing equation of the problem while matching discrete surface measurements. A crucial eikonal regularization term in the loss function ensures that the material density transitions between 0 and 1 over a prescribed length scale and avoids settling on intermediate values. By the end of the optimization, the converged material distribution reveals the location and shapes of the hidden structures.
to vanish on a predefined boundary. Rather, the zero level-set of \(\phi\) evolves freely during the training process in such a way that the corresponding material distribution \(\rho=\text{sigmoid}(\phi/\delta)\) and physical quantities \(\mathbf{\psi}\) minimize the total loss (3), eventually revealing the material boundaries delineating the hidden voids or inclusions.

### Setup of numerical experiments

We evaluate our TO framework on a range of challenging test cases involving different numbers and shapes of hidden structures and various materials (Methods, Tabs. I and II). As a substitute for real experiments, we use the finite-element method (FEM) software Abaqus to compute the deformed shape of the boundary of the elastic structure and generate the measurement data for each case (Methods, section "FEM simulations"). Using this measurement data, we run our TO framework to discover the number, locations, and shapes of the hidden voids or rigid inclusions (for implementation and training details, see Methods, section "Architecture and training details"). We then compare the obtained results with the ground truth -- the voids or inclusions originally fed into Abaqus -- to assess the efficacy of our framework.

### Elastic matrix experiments

We first apply our framework to cases involving a linear elastic matrix (Fig. 1a) containing voids (cases 1, 3, 8, 10, 15, 17 in Methods, Tab. 1). As the various loss components are minimized during training (Fig. 4a), the material density \(\rho\) evolves and splits in a way that progressively reveals the number, locations, and shapes of the hidden voids (Extended Data Fig. 1 and Supplementary Movie 1), without advance knowledge of their topology. By the end of the training, the transition regions where the material density goes from 0 to 1 have uniform thickness along all internal boundaries (Fig. 4d), thanks to the eikonal regularization that encourages the level-set gradient \(\nabla\phi\) to have unit norm in a band along the material boundaries \(\phi=0\) (Fig. 4b,c). The agreement between the final inferred shapes and the ground truth is remarkable, with our framework able to recover intricate details such as the three lobes and the concave surfaces of the star-shaped void (Fig. 4d, second from left), or the exact aspect ratio and location of a thin slit (Fig. 4d, third from left). The stress and strain fields of the deformed matrix are also obtained as a byproduct of the solution process (Extended Data Fig. 2). The only case that is not completely identified is the U-shaped void (Fig. 4d, first from right), a result of the miniscule influence of the inner lobe on the outer surface displacements due to its low level of strain and stress (Extended Data Fig. 8). Finally, our framework maintains accurate results when reducing the number of surface measurement points or restricting measurements to a few surfaces (Extended Data Figs. 5-7).

Next, we consider cases involving linear elastic and rigid inclusions in the linear elastic matrix (cases 4, 5, 6, 11, 12, 13 in Methods, Tab. 1). Our framework successfully identifies the inclusions in almost all cases (Fig. 4e). Inferred displacements and stresses of the deformed matrix (Extended Data Fig. 3) confirm the intuition that voids or soft inclusions soften the matrix while stiff or

Figure 3: **Eikonal regularization of the material density.****a**, A random level-set function \(\phi\) yields a material density \(\rho=\text{sigmoid}(\phi/\delta)\) with large regions of values between 0 and 1, due to the nonuniformity of the gradient \(|\nabla\phi|\) along the material boundaries defined by the zero level-set of \(\phi\) (black lines). **b**, Constraining \(\phi\) to solve the eikonal equation \(\nabla\phi=1\) in a narrow band \(\Omega_{\text{eik}}\) of thickness \(w\) (edges depicted by dashed lines) along the material boundaries results in a uniform transition thickness of \(\rho\) from 0 to 1, without large regions of intermediate density values. **c**, The loss \(\mathcal{L}_{\text{eik}}\) implements the eikonal regularization in the PINN-based TO framework by penalizing deviations away from the constraint \(|\nabla\phi|=1\) on a subset of collocation points \(\Omega_{\text{eik}}^{d}\subset\Omega_{d}\) that approximates the true narrow band \(\Omega_{\text{eik}}\).
Figure 4: **Identification of voids and inclusions in elastic matrices.** A linear elastic matrix containing voids (**a-d**): **a**, The various loss components that enforce the solution to match the surface measurement data, satisfy the governing equations, and obey the eikonal regularization, are being minimized during the training process. **b,c**, The final level-set function \(\phi\) and its gradient magnitude \(|\nabla\phi|\) show the effect of the eikonal regularization, making \(\phi\) a signed distance function in narrow band along the interface. **d**, The final material density \(\rho\) reveals the number, locations, and shapes of the hidden voids, which are compared with the ground truth shown in dotted white lines. **e**, The final material density predictions in the case of a linear elastic matrix containing soft, stiff or rigid inclusions. **f**, The final material density predictions in the case of a nonlinear hyperelastic matrix containing voids subject to large stretches.
rigid inclusions harden the matrix. The U-shaped soft and stiff elastic inclusions (Fig. 4e, second and third from right) are better detected by the framework than their void or rigid counterparts (Fig. 4d, first from right and Fig. 4e, first from right), since an elastic inclusion induces some strain and stress on the inner lobe (Extended Data Fig. 8).

Finally, we consider cases involving a soft, incompressible Neo-Hookean hyperelastic matrix with the same void shapes considered previously (cases 2, 7, 9, 14, 16, 18 in Methods, Tab. 1). The geometries are identified equally well (Fig. 4f) in this large deformation regime (Extended Data Fig. 4) as with linear elastic materials, which illustrates the ability of the framework to cope with nonlinear governing equations without any added complexity in the formulation or the implementation.

### Elastic layer experiments

We finally apply our framework to the periodic elastic layer (Fig. 1b), where a linear elastic material covers a hidden rigid substrate (cases 20, 21, 22 in Methods, Tab. 2).

Contrary to the matrix problem, this setup only provides access to measurements on the top surface, and the hidden geometry to be discovered is not completely surrounded by the elastic material. Our TO framework is nevertheless able to detect the correct depths and shapes of the hidden substrates (Fig. 5). This example demonstrates the versatility of the framework in adapting to various problem setups.

## Discussion

As with any TO method relying on a material density field to parameterize the geometry, the success of our PINN-based framework hinges on the presence of an appropriate regularization mechanism to penalize intermediate density values. Although we have shown that our novel eikonal regularization leads to consistently accurate results, other regularization approaches have been employed in classical adjoint-based TO methods [51, 18]. These include the total variation dimishing (TVD) regularization [39, 10] that penalizes the \(L_{1}\) norm of the density gradient \(\nabla\rho\), the explicit penalization regularization [3] that penalizes the integral over the domain of \(\rho(1-\rho)\), and the Solid Isotropic Material with Penalization (SIMP) approach [7] that relates material properties such as the shear modulus and the material density through a power-law with exponent \(p\). The latter is the most popular regularization mechanism in structural optimization [8]. However, when implemented in our PINN-based framework for the detection of hidden geometries, these methods yield inferior results to the eikonal regularization (Fig. 6). Indeed, we compare all four approaches on a challenging test case involving a linear elastic rectangular matrix pulled from the top and bottom and containing soft inclusions in the shape of the letters M, I, and T (case 19, Tab. 1). The measurements consist of the displacement along the outer boundary, similar to the previous square matrix examples. We consider different values of the regularization weight \(\lambda_{\text{reg}}\) (for the eikonal, TVD and explicit penalization regularizations) and the exponent \(p\) (for the SIMP regularization), and solve the inverse problem using four random initializations of the neural networks in each case. Not only was the eikonal regularization the only one to find the right shapes, it did so over three orders of magnitude of \(\lambda_{\text{reg}}\), demonstrating a desirable robustness with respect to \(\lambda_{\text{reg}}\) (Fig. 6).

Thanks to the flexibility of the PINN framework, adapting our approach to various scenarios involving partial information, three-dimensional geometries, or other types of noninvasive imaging experiments (e.g. using thermal [6], acoustic [14], electric [12], or magnetic [37] loading) should be straightforward. As an illustration, we identify a hidden inclusion in a nonlinearly conducting matrix using partially unknown thermal loading (Supplementary Information), mimicking an inaccessible surface. This experiment reveals our method's ability to generate good results without further modifications even when the forward problem is ill-posed, which would require including the unknown boundary condition as an additional op

Figure 5: **Identification of substrate shape underneath a periodic linear elastic layer.****a**, The various loss components that enforce the solution to match the surface measurement data, satisfy the governing equations, and obey the eikonal regularization, are being minimized during the training process. **b,c**, The final level-set function \(\phi\) and its gradient magnitude \(|\nabla\phi|\) show the effect of the eikonal regularization, which makes \(\phi\) a signed distance function in narrow band along the material boundary. **d**, The final material density \(\rho\) reveals the shape of the buried rigid substrate.
timization variable in a classical adjoint-based approach.

In conclusion, we have presented a PINN-based TO framework with a novel eikonal regularization, which we have applied to the noninvasive detection of hidden inclusions. By representing the geometry through a material density field combined with a novel eikonal regularization, our framework is able to discover the number, shapes and locations of hidden structures, without any prior knowledge required regarding the number or the types of shapes to expect. Finally, the idea of parameterizing geometries of arbitrary topologies with a material density field regularized with the eikonal constraint opens a pathway for PINNs to be applied to a wide range of design optimization problems constrained by physical governing equations. These include, for instance, the design of lenses that achieve targeted optical properties [38, 41] or the design of structures and metamaterials that exhibit desirable mechanical, acoustic, or thermal properties [32, 8, 30].

Fig. 6: **Comparison between eikonal regularization and alternative regularizations.****a**, The eikonal regularization achieves a high IoU (intersection over union, a geometry detection accuracy metric equal to 1 in the perfect case) above 0.95 for any value of the regularization weight \(\lambda_{\text{reg}}\) within a range spanning three orders of magnitude. The results are consistent over 4 random initializations of the neural networks parameters, with the circles reporting the average value and the shade reporting the highest and lowest values. By contrast, the total variation diminishing (TVD), explicit penalization, and Solid Isotropic Material with Penalization (SIMP) regularizations never exceed IoU values above 0.93, with larger variability among realizations. **b**, The final material density \(\rho\) obtained with each regularization mechanism for various values of the regularization weight \(\lambda_{\text{reg}}\) or exponent \(p\) (shown in **a** by the shaded areas) demonstrates the efficacy of the eikonal regularization. The evolution of the solution during training using the eikonal regularization and \(\lambda_{\text{reg}}=1\) is shown in Supplementary Movie 2.
## Methods

### Governing equations

The two plane-strain elasticity inverse problems considered in this study (Fig. 1) are defined in a two-dimensional domain \(\Omega\subset\mathbb{R}^{2}\) formed by the union of the elastic body \(\mathcal{B}\) and the voids or inclusions. Denoting with \(\mathbf{x}=(x_{1},x_{2})\in\Omega\) the planar spatial coordinates, the hidden geometrical layout of voids or inclusions is characterized by a material density \(\rho(\mathbf{x})\) equal to \(1\) in the body \(\mathcal{B}\) and \(0\) in the voids or inclusions.

#### Small-deformation linear elasticity

We first consider the case where the elastic body and inclusions consist of linear elastic materials, with Young's modulus \(E\) and Poisson's ratio \(\nu\) for the body, and Young's modulus \(\bar{E}\) and Poisson's ratio \(\bar{\nu}\) for the inclusions. Voids and rigid inclusions correspond to the limits \(\bar{E}\to 0\) and \(\bar{E}\to\infty\), respectively. The deformation of the elastic body containing the inclusions is described by a vector field \(\mathbf{\psi}(\mathbf{x})=(\mathbf{u}(\mathbf{x}),\mathbf{\sigma}(\mathbf{x}))\), where \(\mathbf{u}(\mathbf{x})\) is a planar displacement field with components \(u_{i}(\mathbf{x})\) and \(\mathbf{\sigma}(\mathbf{x})\) is a Cauchy stress tensor with components \(\sigma_{ij}(\mathbf{x})\). Indices \(i\) and \(j\) will hereafter always range from \(1\) to \(2\).

The governing PDEs comprise the equilibrium equations

\[\sum_{j}\frac{\partial\sigma_{ij}}{\partial x_{j}}=0,\quad\mathbf{x}\in\Omega, \tag{7}\]

as well as a linear elastic constitutive law \(F(\mathbf{\sigma},\nabla\mathbf{u},\rho)=0\) that we will express in two different but equivalent ways, depending on whether the inclusions are softer or stiffer than the matrix. For voids and soft inclusions, we consider the constitutive law in stress-strain form,

\[\mathbf{\sigma}= \,\rho\left[\lambda\,\mathrm{tr}(\mathbf{\epsilon})\,\mathbf{I}+2\mu \,\mathbf{\epsilon}\right],\] \[(1-\rho)\left[\bar{\lambda}\,\mathrm{tr}(\mathbf{\epsilon})\,\mathbf{ I}+2\bar{\mu}\,\mathbf{\epsilon}\right],\quad\mathbf{x}\in\Omega, \tag{8}\]

where \(\mathbf{\epsilon}=(\nabla\mathbf{u}+\nabla\mathbf{u}^{T})/2\) is the infinitesimal strain tensor, \(\mathrm{tr}(\mathbf{\epsilon})\) denotes its trace, \(\lambda=E\nu/(1+\nu)(1-2\nu)\) and \(\mu=E/2(1+\nu)\) are the Lame constants of the body, and \(\bar{\lambda}=\bar{E}\bar{\nu}/(1+\bar{\nu})(1-2\bar{\nu})\) and \(\bar{\mu}=\bar{E}/2(1+\bar{\nu})\) are the Lame constants of the inclusions. Notice that the case of voids, the stress vanishes in the \(\rho=0\) regions. For stiff and rigid inclusions, we consider the constitutive law in the inverted strain-stress form

\[\mathbf{\epsilon}= \,\rho\left[\frac{1+\nu}{E}\mathbf{\sigma}-\frac{\nu(1+\nu)}{E}\, \mathrm{tr}(\mathbf{\sigma})\,\mathbf{I}\right]\] \[(1-\rho)\left[\frac{1+\bar{\nu}}{\bar{E}}\mathbf{\sigma}-\frac{\bar{ \nu}(1+\bar{\nu})}{\bar{E}}\,\mathrm{tr}(\mathbf{\sigma})\,\mathbf{I}\right],\quad \mathbf{x}\in\Omega, \tag{9}\]

where \(\mathrm{tr}(\mathbf{\sigma})\) is the trace of the stress tensor. This relation differs from the three-dimensional one due to the plane strain assumption. Notice that the case of rigid inclusions, the strain vanishes in the \(\rho=0\) regions.

The boundary conditions on \(\partial\Omega\) and surface displacement measurement locations \(\partial\Omega^{m}\) are different in the two problems. For the elastic matrix (Fig. 1a), the domain is \(\Omega=[-0.5,0.5]\times[-0.5,0.5]\) and the boundary conditions are

\[\mathbf{\sigma}(\mathbf{x})\mathbf{n}(\mathbf{x}) =-P_{o}\mathbf{e}_{1}, \mathbf{x}\in\{-0.5,0.5\}\times[-0.5,0.5], \tag{10a}\] \[\mathbf{\sigma}(\mathbf{x})\mathbf{n}(\mathbf{x}) =\mathbf{0}, \mathbf{x}\in[-0.5,0.5]\times\{-0.5,0.5\}. \tag{10b}\]

The measurement locations \(\partial\Omega^{m}\) are distributed along the entire external boundary \(\partial\Omega\). In the case of the M, I, T inclusions (case 19, Tab. 1), the boundary conditions (10) are changed to account for the fact that the matrix is pulled from the top and bottom boundaries and covers the domain \(\Omega=[-1,1]\times[-0.5,0.5]\). For the elastic layer (Fig. 1b), the domain is \(\Omega=[0,1]\times[-0.5,0]\) and the boundary conditions are

\[\mathbf{\sigma}(\mathbf{x})\mathbf{n}(\mathbf{x}) =-P_{o}\mathbf{e}_{2}, \mathbf{x}\in[0,1]\times\{0\}, \tag{11a}\] \[\mathbf{u} =\mathbf{0}, \mathbf{x}\in[0,1]\times\{-0.5\}, \tag{11b}\]

as well as periodic for the displacement and traction on \(\mathbf{x}\in\{0,1\}\times[-0.5,0]\). The measurement locations \(\partial\Omega^{m}\) are distributed along the top surface \(\partial\Omega_{t}=[0,1]\times\{0\}\).

The geometry identification problem that we solve can then be stated as follows. Given surface displacement measurements \(\mathbf{u}_{i}^{m}\) at locations \(\mathbf{x}_{i}\in\partial\Omega^{m}\), find the distribution of material density \(\rho\) in \(\Omega\) such that the difference between the predicted and measured surface displacements vanish, that is,

\[\mathbf{u}(\mathbf{x}_{i})=\mathbf{u}_{i}^{m},\quad\mathbf{x}_{i}\in\partial \Omega^{m}. \tag{12}\]

The predicted displacement field must satisfy the equilibrium equation (7), the constitutive relation (8) or (9), and the boundary conditions (10) or (11).

#### Large-deformation nonlinear hyperelasticity

Next, we consider the case where the elastic body consists of an incompressible Neo-Hookean hyperelastic material with shear modulus \(\mu\). We now have to distinguish between the reference (undeformed) and current (deformed) configurations. We denote by \(\mathbf{x}=(x_{1},x_{2})\in\Omega\) and \(\mathbf{y}=(y_{1},y_{2})\in\Omega^{*}\) the coordinates in the reference and deformed configurations, respectively, with \(\Omega^{*}\) the deformed image of \(\Omega\). The displacement field \(\mathbf{u}(\mathbf{x})\) with components \(u_{i}(\mathbf{x})\) moves an initial position \(\mathbf{x}\in\Omega\) into its current location \(\mathbf{y}=\mathbf{x}+\mathbf{u}(\mathbf{x})\in\Omega^{*}\). In order to formulate the governing equations and boundary conditions in the reference configuration \(\Omega\), we need to introduce the first Piola-Kirchhoff stress tensor \(\mathbf{S}(\mathbf{x})\) with components \(S_{ij}(\mathbf{x})\). Unlike the Cauchy stress tensor, the first Piola-Kirchhoff stress tensor is defined in \(\Omega\) and is not
symmetric. The deformation of the elastic body is then described by the vector field \(\mathbf{\psi}(\mathbf{x})=(\mathbf{u}(\mathbf{x}),\mathbf{S}(\mathbf{x}),p(\mathbf{x}))\) defined over \(\Omega\), where \(p(\mathbf{x})\) is a pressure field that serves to enforce the incompressibility constraint.

The equilibrium equations are

\[\sum_{j}\frac{\partial S_{ij}}{\partial x_{j}}=0,\quad\mathbf{x}\in\Omega, \tag{13}\]

where the derivatives in \(\nabla_{\mathbf{x}}\) are taken with respect to the reference coordinates \(\mathbf{x}\). We only consider the presence of voids so that the nonlinear constitutive law \(F(\mathbf{S},\nabla_{\mathbf{x}}\mathbf{u},p,\rho)=0\) is simply expressed as

\[\mathbf{S}=\rho\left[-p\mathbf{F}^{-T}+\mu\mathbf{F}\right],\quad\mathbf{x}\in\Omega, \tag{14}\]

where \(\mathbf{F}(\mathbf{x})=\mathbf{I}+\nabla_{\mathbf{x}}\mathbf{u}(\mathbf{x})\) is the deformation gradient tensor. Notice that the stress vanishes in the \(\rho=0\) regions. Finally, we have the incompressibility constraint

\[\rho\left[\det(\mathbf{F})-1\right]=0,\quad\mathbf{x}\in\Omega, \tag{15}\]

which turns itself off in the \(\rho=0\) regions since voids do not deform in a way that preserves volume.

We only treat the matrix problem (Fig. 1a) in this hyperelastic case. The domain is \(\Omega=[-0.5,0.5]\times[-0.5,0.5]\) and the boundary conditions are

\[\mathbf{S}(\mathbf{x})\mathbf{n}_{0}(\mathbf{x}) =-P_{o}\mathbf{e}_{1}, \mathbf{x}\in\{-0.5,0.5\}\times[-0.5,0.5], \tag{16a}\] \[\mathbf{S}(\mathbf{x})\mathbf{n}_{0}(\mathbf{x}) =\mathbf{0}, \mathbf{x}\in[-0.5,0.5]\times\{-0.5,0.5\}. \tag{16b}\]

As in the linear elastic case, the measurement locations \(\partial\Omega^{m}\) are distributed along the entire external boundary \(\partial\Omega\).

The geometry identification problem can then be stated identically as in the linear elastic case. This time, the predicted displacement field must satisfy the equilibrium equation (13), the constitutive relation (14), the incompressibility condition (15), and the boundary conditions (16).

#### iv.2.2 Rescaling

The various physical quantities involved in the elasticity inverse problem span a wide range of scales; for instance, displacements may be orders of magnitude smaller than the length scale associated with the geometry. Thus, we rescale all physical quantities into nondimensional values of order one, as also done in [29]. Lengths are rescaled with the width \(L\) of the elastic matrix or elastic layer, tractions and stresses with the magnitude \(P_{o}\) of the applied traction at the boundaries, and displacements with the ratio \(LP_{o}/E\), where \(E\) is the Young's modulus of the elastic material (in the hyperelastic case, we use the equivalent Young's modulus \(E=3\mu\), where \(\mu\) is the shear modulus of the hyperelastic material). This rescaling is critical to enable the neural networks underlying our framework to handle elasticity problems across a wide range of material moduli and applied loads.

### Solution methodology

Here, we describe in detail the application of our TO framework to the solution of the two plane-strain elasticity inverse problems formulated in the introduction. We will treat separately the small-deformation linear elasticity case and the large-deformation hyperelasticity case.

#### iv.2.3 Small-deformation linear elasticity

Since the problem is described by the physical quantities \(\mathbf{\psi}=(u_{1},u_{2},\sigma_{11},\sigma_{22},\sigma_{12})\), we introduce the neural network approximations

\[u_{1}(\mathbf{x}) =\bar{u}_{1}(\mathbf{x};\mathbf{\theta}_{1}), \tag{17a}\] \[u_{2}(\mathbf{x}) =\bar{u}_{2}(\mathbf{x};\mathbf{\theta}_{2}),\] (17b) \[\sigma_{11}(\mathbf{x}) =\bar{\sigma}_{11}(\mathbf{x};\mathbf{\theta}_{3}),\] (17c) \[\sigma_{22}(\mathbf{x}) =\bar{\sigma}_{22}(\mathbf{x};\mathbf{\theta}_{4}),\] (17d) \[\sigma_{12}(\mathbf{x}) =\bar{\sigma}_{12}(\mathbf{x};\mathbf{\theta}_{5}),\] (17e) \[\phi(\mathbf{x}) =\bar{\phi}(\mathbf{x};\mathbf{\theta}_{\phi}). \tag{17f}\]

The last equation represents the level-set neural network, which defines the material density as \(\rho(\mathbf{x})=\bar{\rho}(\mathbf{x};\mathbf{\theta}_{\phi})=\text{sigmoid}( \bar{\phi}(\mathbf{x};\mathbf{\theta}_{\phi})/\delta)\). We then formulate the loss function (3) by specializing the loss term expressions presented in results section to the linear elasticity problem. Omitting the \(\mathbf{\theta}\)'s for notational simplicity, we obtain

\[\mathcal{L}_{\text{meas}}(\mathbf{\theta}_{\mathbf{\psi}}) =\frac{1}{|\partial\Omega^{m}|}\sum_{\mathbf{x}_{i}\in\partial \Omega^{m}}|\bar{\mathbf{u}}(\mathbf{x}_{i})-\mathbf{u}_{i}^{m}|^{2}, \tag{18a}\] \[\mathcal{L}_{\text{gov}}(\mathbf{\theta}_{\mathbf{\psi}},\mathbf{\theta}_{ \phi}) =\frac{1}{|\Omega^{d}|}\sum_{\mathbf{x}_{i}\in\Omega^{d}}|\mathbf{r }_{\text{eq}}(\bar{\mathbf{\sigma}}(\mathbf{x}_{i}))|^{2}\] \[\quad+\frac{1}{|\Omega^{d}|}\sum_{\mathbf{x}_{i}\in\Omega^{d}}| \mathbf{r}_{\text{cr}}(\bar{\mathbf{u}}(\mathbf{x}_{i}),\bar{\mathbf{\sigma}}( \mathbf{x}_{i}),\bar{\rho}(\mathbf{x}_{i}))|^{2},\] (18b) \[\mathcal{L}_{\text{eik}}(\mathbf{\theta}_{\phi}) =\frac{1}{|\Omega^{d}_{\text{eik}}|}\sum_{\mathbf{x}_{i}\in \Omega^{d}_{\text{eik}}}\left(|\nabla\bar{\phi}(\mathbf{x}_{i})|-1\right)^{2}, \tag{18c}\]

where \(\bar{\mathbf{u}}=(\bar{u}_{1},\bar{u}_{2})\) and \(\bar{\mathbf{\sigma}}\) has components \(\bar{\sigma}_{i,j}\), \(i,j=1,2\). In (18b), the terms \(\mathbf{r}_{\text{eq}}\) and \(\mathbf{r}_{\text{cr}}\) refer to the residuals of the equilibrium equation (7) and the constitutive relation (8) or (9). The eikonal loss term is problem-independent and therefore identical to (6).

We note that instead of defining neural network approximations for the displacements and the stresses, we could define neural network approximations solely for the displacements, that is, \(\mathbf{\psi}=(u_{1},u_{2})\). In this case, the loss term (18b) would only include the residual of the equilibrium equation (7), in which the stress components would
be directly expressed in terms of the displacements and the material distribution using the constitutive relation (8). However, several recent studies [21, 25, 26, 29, 47, 48] have shown that the mixed formulation adopted in the present work results in superior accuracy and training performance, which could partly be explained by the fact that only first-order derivatives of the neural network outputs are involved since the displacements and stresses are only differentiated to first oder in (7) and (8). In our case, the mixed formulation holds the additional advantage that it enables us to treat stiff and rigid inclusions using the inverted constitutive relation (9) instead of (8). Finally, the mixed formulation allows us to directly integrate both displacement and traction boundary conditions into the output of the neural network approximations, as we describe in the next paragraph.

We design the architecture of the neural networks in such a way that they inherently satisfy the boundary conditions, treating the latter as hard constraints [17, 53]. For the elastic matrix, we do this through the transformations

\[\bar{u}_{1}(\mathbf{x};\mathbf{\theta}_{1}) =\bar{u}^{\prime}_{1}(\mathbf{x};\mathbf{\theta}_{1}), \tag{19a}\] \[\bar{u}_{2}(\mathbf{x};\mathbf{\theta}_{2}) =\bar{u}^{\prime}_{2}(\mathbf{x};\mathbf{\theta}_{2}),\] (19b) \[\bar{\sigma}_{11}(\mathbf{x};\mathbf{\theta}_{3}) =(x-0.5)(x+0.5)\,\bar{\sigma}^{\prime}_{11}(\mathbf{x};\mathbf{\theta }_{3})+P_{o},\] (19c) \[\bar{\sigma}_{22}(\mathbf{x};\mathbf{\theta}_{4}) =(y-0.5)(y+0.5)\,\bar{\sigma}^{\prime}_{22}(\mathbf{x};\mathbf{\theta }_{4}),\] (19d) \[\bar{\sigma}_{12}(\mathbf{x};\mathbf{\theta}_{5}) =(x-0.5)(x+0.5)\cdot\] \[(y-0.5)(y+0.5)\,\bar{\sigma}^{\prime}_{12}(\mathbf{x};\mathbf{\theta }_{5}),\] (19e) \[\bar{\phi}(\mathbf{x};\mathbf{\theta}_{\phi}) =(x-0.5)(x+0.5)\cdot\] \[(y-0.5)(y+0.5)\,\bar{\phi}^{\prime}(\mathbf{x};\mathbf{\theta}_{\phi} )+w, \tag{19f}\]

where the quantities with a prime denote the raw output of the neural network. In this way, the neural network approximations defined in (17) obey by construction the boundary conditions (10). Further, since we know that the elastic material is present all along the outer surface \(\partial\Omega\), we define \(\bar{\phi}\) so that \(\phi=w\) on \(\partial\Omega\), which ensures that \(\rho=\text{sigmoid}(\phi/\delta)\simeq 1\) on \(\partial\Omega\) (recall that \(w\) is such that \(\text{sigmoid}(w/2\delta)\simeq 1\)). In the case of the M, I, T inclusions, these transformations are changed to reflect the fact that the matrix is wider and pulled from the top and bottom. For the periodic elastic layer, we introduce the transformations

\[\bar{u}_{1}(\mathbf{x};\mathbf{\theta}_{1}) =(y+0.5)\,\bar{u}^{\prime}_{1}(\cos x,\sin x,y;\mathbf{\theta}_{1}), \tag{20a}\] \[\bar{u}_{2}(\mathbf{x};\mathbf{\theta}_{2}) =(y+0.5)\,\bar{u}^{\prime}_{2}(\cos x,\sin x,y;\mathbf{\theta}_{2}),\] (20b) \[\bar{\sigma}_{11}(\mathbf{x};\mathbf{\theta}_{3}) =\bar{\sigma}^{\prime}_{11}(\cos x,\sin x,y;\mathbf{\theta}_{3}),\] (20c) \[\bar{\sigma}_{22}(\mathbf{x};\mathbf{\theta}_{4}) =y\,\bar{\sigma}^{\prime}_{22}(\cos x,\sin x,y;\mathbf{\theta}_{4})-P _{o},\] (20d) \[\bar{\sigma}_{12}(\mathbf{x};\mathbf{\theta}_{5}) =y\,\bar{\sigma}^{\prime}_{12}(\cos x,\sin x,y;\mathbf{\theta}_{5}),\] (20e) \[\bar{\phi}(\mathbf{x};\mathbf{\theta}_{\phi}) =y(y+0.5)\,\bar{\phi}^{\prime}(\cos x,\sin x,y;\mathbf{\theta}_{\phi})\] \[\quad+w(4y+1), \tag{20f}\]

so that the neural network approximations defined in (17) obey by construction the boundary conditions (11) and are periodic along the \(x\) direction. Further, since we know that the elastic material is present all along the top surface \(y=0\) and the rigid substrate is present all along the bottom surface \(y=-0.5\), we define \(\bar{\phi}\) so that \(\phi=w\) for \(y=0\) and \(\phi=-w\) for \(y=-0.5\), which ensures that \(\rho=\text{sigmoid}(\phi/\delta)\simeq 1\) for \(y=0\) and \(\rho\simeq 0\) for \(y=-0.5\).

#### Large-deformation hyperelasticity

The problem is now described by the physical quantities \(\mathbf{\psi}=(u_{1},u_{2},S_{11},S_{22},S_{12},S_{21},p)\). We therefore introduce the neural network approximations

\[u_{1}(\mathbf{x}) =\bar{u}_{1}(\mathbf{x};\mathbf{\theta}_{1}), \tag{21a}\] \[u_{2}(\mathbf{x}) =\bar{u}_{2}(\mathbf{x};\mathbf{\theta}_{2}),\] (21b) \[S_{11}(\mathbf{x}) =\bar{S}_{11}(\mathbf{x};\mathbf{\theta}_{3}),\] (21c) \[S_{22}(\mathbf{x}) =\bar{S}_{22}(\mathbf{x};\mathbf{\theta}_{4}),\] (21d) \[S_{12}(\mathbf{x}) =\bar{S}_{12}(\mathbf{x};\mathbf{\theta}_{5}),\] (21e) \[S_{21}(\mathbf{x}) =\bar{S}_{21}(\mathbf{x};\mathbf{\theta}_{6}),\] (21f) \[p(\mathbf{x}) =\bar{p}(\mathbf{x};\mathbf{\theta}_{7}),\] (21g) \[\phi(\mathbf{x}) =\bar{\phi}(\mathbf{x};\mathbf{\theta}_{\phi}), \tag{21h}\]

and the material distribution is given by \(\rho(\mathbf{x})=\bar{p}(\mathbf{x};\mathbf{\theta}_{\phi})=\text{sigmoid}(\bar{ \phi}(\mathbf{x};\mathbf{\theta}_{\phi})/\delta)\). We then formulate the loss function (3) by specializing the loss term expressions presented in Section 4 to the linear elasticity problem, using the governing equations given in Appendix. Omitting the \(\mathbf{\theta}\)'s for notational simplicity, we obtain

\[\mathcal{L}_{\text{meas}}(\mathbf{\theta}_{\mathbf{\psi}}) =\frac{1}{|\partial\Omega^{m}|}\sum_{\mathbf{x}_{i}\in\partial \Omega^{m}}|\bar{\mathbf{u}}(\mathbf{x}_{i})-\mathbf{u}_{i}^{m}|^{2}, \tag{22a}\] \[\mathcal{L}_{\text{gov}}(\mathbf{\theta}_{\mathbf{\psi}},\mathbf{\theta}_{ \phi}) =\frac{1}{|\Omega^{d}|}\sum_{\mathbf{x}_{i}\in\Omega^{d}}|\mathbf{r}_{ \text{eq}}(\bar{\mathbf{S}}(\mathbf{x}_{i}))|^{2}\] \[\quad+\frac{1}{|\Omega^{d}|}\sum_{\mathbf{x}_{i}\in\Omega^{d}}| \mathbf{r}_{\text{cr}}(\bar{\mathbf{u}}(\mathbf{x}_{i}),\bar{\mathbf{S}}( \mathbf{x}_{i}),\bar{p}(\mathbf{x}_{i}),\bar{\rho}(\mathbf{x}_{i}))|^{2}\] \[\quad+\frac{1}{|\Omega^{d}|}\sum_{\mathbf{x}_{i}\in\Omega^{d}}| \mathbf{r}_{\text{inc}}(\bar{\mathbf{u}}(\mathbf{x}_{i}),\bar{\rho}(\mathbf{x}_ {i}))|^{2},\] (22b) \[\mathcal{L}_{\text{eik}}(\mathbf{\theta}_{\phi}) =\frac{1}{|\Omega^{d}_{\text{eik}}|}\sum_{\mathbf{x}_{i}\in \Omega^{d}_{\text{eik}}}\left(|\nabla\bar{\phi}(\mathbf{x}_{i})|-1\right)^{2}, \tag{22c}\]

where \(\bar{\mathbf{u}}=(\bar{u}_{1},\bar{u}_{2})\) and \(\bar{\mathbf{S}}\) has components \(\bar{S}_{i,j}\), \(i,j=1,2\). In (22b), the terms \(\mathbf{r}_{\text{eq}}\), \(\mathbf{r}_{\text{cr}}\), and \(\mathbf{r}_{\text{inc}}\) refer to the residuals of the equilibrium equation (13), the constitutive relation (14), and the incompressibility constraint (15). The eikonal loss term is problem-independent and therefore identical to (6).

As in the linear elasticity case, we design the architecture of the neural networks in such a way that they inherently satisfy the boundary conditions. For the elastic matrix problem,

\[\bar{u}_{1}(\mathbf{x};\mathbf{\theta}_{1}) =\bar{u}^{\prime}_{1}(\mathbf{x};\mathbf{\theta}_{1}), \tag{23a}\] \[\bar{u}_{2}(\mathbf{x};\mathbf{\theta}_{2}) =\bar{u}^{\prime}_{2}(\mathbf{x};\mathbf{\theta}_{2}), \tag{23b}\]
\[\bar{S}_{11}(\mathbf{x};\mathbf{\theta}_{3}) =(x-0.5)(x+0.5)\,\bar{S}^{\prime}_{11}(\mathbf{x};\mathbf{\theta}_{3})+P _{o}, \tag{23c}\] \[\bar{S}_{22}(\mathbf{x};\mathbf{\theta}_{4}) =(y-0.5)(y+0.5)\,\bar{S}^{\prime}_{22}(\mathbf{x};\mathbf{\theta}_{4}),\] (23d) \[\bar{S}_{12}(\mathbf{x};\mathbf{\theta}_{5}) =(y-0.5)(y+0.5)\,\bar{S}^{\prime}_{12}(\mathbf{x};\mathbf{\theta}_{5}),\] (23e) \[\bar{S}_{21}(\mathbf{x};\mathbf{\theta}_{6}) =(x-0.5)(x+0.5)\,\bar{S}^{\prime}_{21}(\mathbf{x};\mathbf{\theta}_{6}),\] (23f) \[\bar{p}(\mathbf{x};\mathbf{\theta}_{7}) =\bar{p}^{\prime}(\mathbf{x};\mathbf{\theta}_{7}),\] (23g) \[\bar{\phi}(\mathbf{x};\mathbf{\theta}_{\phi}) =(x-0.5)(x+0.5)\cdot\] \[\quad(y-0.5)(y+0.5)\,\bar{\phi}^{\prime}(\mathbf{x};\mathbf{\theta}_{ \phi})+w, \tag{23h}\]

where the quantities with a prime denote the raw output of the neural network. In this way, the neural network approximations defined in (21) obey by construction the boundary conditions of the problem. As before, since we know that the elastic material is present all along the outer surface \(\partial\Omega\), we define \(\bar{\phi}\) so that \(\phi=w\) on \(\partial\Omega\), which ensures that \(\rho=\mathrm{sigmoid}(\phi/\delta)\simeq 1\) on \(\partial\Omega\).

### Architecture and training details

Here, we provide implementation details regarding the architecture of the deep neural networks, the training procedure and corresponding parameter values.

#### Neural network architecture

State variable fields of the form \(\psi(\mathbf{x})\) are approximated using deep fully-connected neural networks that map the location \(\mathbf{x}\) to the corresponding value of \(\psi\) at that location. This map can be expressed as \(\psi(\mathbf{x})=\bar{\psi}(\mathbf{x};\mathbf{\theta})\), and is defined by the sequence of operations

\[\mathbf{z}^{0} =\mathbf{x}, \tag{24a}\] \[\mathbf{z}^{k} =\sigma(\mathbf{W}^{k}\mathbf{z}^{k-1}+\mathbf{b}^{k}),\quad 1 \leq k\leq\ell-1,\] (24b) \[\psi =\mathbf{z}^{\ell} =\mathbf{W}^{\ell}\mathbf{z}^{\ell-1}+\mathbf{b}^{\ell}. \tag{24c}\]

The input \(\mathbf{x}\) is propagated through \(\ell\) layers, all of which (except the last) take the form of a linear operation composed with a nonlinear transformation. Each layer outputs a vector \(\mathbf{z}^{k}\in\mathbb{R}^{q_{k}}\), where \(q_{k}\) is the number of 'neurons', and is defined by a weight matrix \(\mathbf{W}^{k}\in\mathbb{R}^{q_{k}\times q_{k-1}}\), a bias vector \(\mathbf{b}^{k}\in\mathbb{R}^{q_{k}}\), and a nonlinear activation function \(\sigma(\cdot)\). Finally, the output of the last layer is assigned to \(\psi\). The weight matrices and bias vectors, which parametrize the map from \(\mathbf{x}\) to \(\psi\), form a set of trainable parameters \(\mathbf{\theta}=\{\mathbf{W}^{k},\mathbf{b}^{k}\}_{k=1}^{\ell}\).

The choice of the nonlinear activation function \(\sigma(\cdot)\) and the initialization procedure for the trainable parameters \(\mathbf{\theta}\) are both important factors in determining the performance of neural networks. While the tanh function has been a popular candidate in the context of PINNs [35], recent works by Refs. [52; 55] have shown that using sinusoidal activation functions can lead to improved training performance by promoting the emergence of small-scale features. In this work, we select the sinusoidal representation network (SIREN) architecture from Ref. [52], which combines the use of the sine as an activation function with a specific way to initialize the trainable parameters \(\mathbf{\theta}\) that ensures that the distribution of the input to each sine activation function remains unchanged over successive layers. Specifically, each component of \(\mathbf{W}^{k}\) is uniformly distributed between \(-\sqrt{6/q_{k}}\) and \(\sqrt{6/q_{k}}\) where \(q_{k}\) is the number of neurons in layer \(k\), and \(\mathbf{b}^{k}=\mathbf{0}\), for \(k=1,\ldots,\ell\). Further, the first layer of the SIREN architecture is \(\mathbf{z}^{1}=\sigma(\omega_{0}\mathbf{W}^{1}\mathbf{z}^{0}+\mathbf{b}^{1})\) instead of (24b), with the extra scalar \(\omega_{0}\) promoting higher-frequency content in the output.

#### Training procedure

We construct the total loss function (3) and train the neural networks in TensorFlow 2. The training is performed using ADAM, a first-order gradient-descent-based algorithm with adaptive step size [31]. In each case, we repeat the training over four random initializations of the neural networks parameters and report the best results. Three tricks resulted in noticeably improved training performance and consistency:

* First, we found that pretraining the level-set neural network \(\phi(\mathbf{x})=\bar{\phi}(\mathbf{x};\mathbf{\theta}_{\phi})\) in a standard supervised setting leads to much more consistent results over different initializations of the neural networks. During this pretraining step, carried out before the main optimization step in which all neural networks are trained to minimize the loss (3), we minimize the mean-square error \[\mathcal{L}_{\mathrm{sup}}(\mathbf{\theta}_{\phi})=\frac{1}{|\Omega^{d}|}\sum_{ \mathbf{x}_{i}\in\Omega^{d}}|\bar{\phi}(\mathbf{x}_{i};\mathbf{\theta}_{\phi})- \phi_{i}|,\] (25) where \(\Omega^{d}\) is the same set of collocation points as in (5), the supervised labels \(\phi_{i}=|\mathbf{x}_{i}|-0.25\) for the elastic matrix, and \(\phi_{i}=y_{i}+0.25\) for the elastic layer. The material density \(\bar{\rho}(\mathbf{x};\mathbf{\theta}_{\phi})=\mathrm{sigmoid}(\bar{\phi}( \mathbf{x};\mathbf{\theta}_{\phi})/\delta)\) obtained at the end of this pretraining step is one outside a circle of radius \(0.25\) centered at the origin for the elastic matrix, and it is one above the horizontal line \(y=-0.25\) for the elastic layer. This choice for the supervised labels is justified by the fact that \(\rho\) is known to be one along the outer boundary of the domain \(\Omega\) for the elastic matrix, and it is known to be one (zero) along the top (bottom) boundary of \(\Omega\) for the elastic layer.
* Second, during the main optimization in which all neural networks are trained to minimize the loss (3), we evaluate the loss component \(\mathcal{L}_{\mathrm{gov}}\) in (5) using a different subset, or mini-batch, of residual points from \(\Omega^{d}\) at every iteration. Such a mini-batching approach has been reported to improve the convergence of the PINN training process [54; 15], corroborating our own observations. In our case, we choose to divide the set \(\Omega^{d}\) into
10 different mini-batches of size \(|\Omega^{d}|/10\), which are then employed sequentially to evaluate \(\mathcal{L}_{\text{gov}}\) during each subsequent gradient update \[\mathbf{\theta}_{\mathbf{\psi}}^{k+1} =\mathbf{\theta}_{\mathbf{\psi}}^{k}-\alpha_{\mathbf{\psi}}(k)\nabla_{\mathbf{\theta }_{\mathbf{\psi}}}\mathcal{L}(\mathbf{\theta}_{\mathbf{\psi}}^{k},\mathbf{\theta}_{\phi}^{k}),\] (26a) \[\mathbf{\theta}_{\phi}^{k+1} =\mathbf{\theta}_{\phi}^{k}-\alpha_{\phi}(k)\nabla_{\mathbf{\theta}_{\phi }}\mathcal{L}(\mathbf{\theta}_{\mathbf{\psi}}^{k},\mathbf{\theta}_{\phi}^{k}).\] (26b) An epoch of training, which is defined as one complete pass through the whole set \(\Omega^{d}\), therefore consists of 10 gradient updates.
* Third, the initial nominal step size \(\alpha_{\mathbf{\psi}}\) governing the learning rate of the physical quantities neural networks is set to be 10 times larger than its counterpart \(\alpha_{\phi}\) governing the learning rate of the level-set neural network. This results in a separation of time scales between the rate of change of the physical quantities neural networks and that of the level-set neural network, which is motivated by the idea that physical quantities should be given time to adapt to a given geometry before the geometry itself changes.

#### Parameter values

The parameter values described below apply to all results presented in this paper.

* **Neural network architecture.** For all cases except the M, I, T inclusions (case 19, Tab. 1), we opted for neural networks with 4 hidden layers of 50 neurons each, which we found to be a good compromise between expressivity and training time. For the M, I, T inclusions, we used 6 hidden layers with 100 neurons each. Further, we choose \(\omega_{0}=10\) as the scalar appearing in the first layer of the SIREN architecture.
* **Collocation and measurements points.** In the square and rectangle elastic matrix problems, we consider that the boundary displacement is measured along each of the four external boundaries at 100 equally-spaced points, which amounts to \(|\partial\Omega^{m}|=400\). In the elastic layer problem, we consider that the boundary displacement is measured along the top boundary at 100 equally-spaced points, which amounts to \(|\partial\Omega^{m}|=100\). For both geometries except the M, I, T inclusions, the set of collocation points \(\Omega^{d}\) consists of 10000 points distributed in \(\Omega\) with a Latin Hypercube Sampling (LHS) strategy, yielding 10 mini-batches containing 1000 points each. For the M, I, T inclusions, \(\Omega^{d}\) consists of 50000 points, yielding 50 mini-batches containing 1000 points each.
* **Training parameters.** The pretraining of the level-set neural network is carried out using the ADAM optimizer with nominal step size \(10^{-3}\) over 800 training epochs, employing the whole set \(\Omega^{d}\) to compute the gradient of \(\mathcal{L}_{\text{sup}}\) at each update step. The main optimization, during which all neural networks are trained to minimize the total loss (3), is carried out using the ADAM optimizer. For the matrix cases except the M, I, T inclusions, we use a total of 150k training epochs starting from a nominal step size \(10^{-4}\) for the level-set neural network and \(10^{-3}\) for the other neural networks. This step size is reduced to \(10^{-4}\) for all neural networks at 60k epochs, and again to \(10^{-5}\) at 120k epochs. The schedule is the same for the elastic layer cases, with the difference that we use a total of 200k training epochs. For the matrix case with the M, I, T inclusions, we use a total of 50k epochs (note that each epoch contains 5 times as many mini-batches as in the other cases) starting from a nominal step size \(10^{-4}\) for the level-set neural network and \(10^{-3}\) for the other neural networks. This step size is reduced to \(10^{-4}\) for all neural networks at 16k epochs, and again to \(10^{-5}\) at 40k epochs. Finally, the scalar weights in the loss (3) are assigned the values \(\lambda_{\text{meas}}=10\), \(\lambda_{\text{gov}}=1\), and \(\lambda_{\text{reg}}=1\) for all cases. We also multiply the second term of \(\mathcal{L}_{\text{gov}}\) in (18b) and (22b) with a scalar weight \(\lambda_{\text{cr}}=10\).

### FEM simulations

The FEM simulations that provide the boundary displacement data and the ground truth are performed in the software Abaqus, using its Standard (implicit) solver. The list of all cases considered in provided in Tab. 1 for the elastic matrix setup (Fig. 1a) and in Tab. 2 for the periodic elastic layer setup (Fig. 1b). Every case is meshed using a linear density of 200 elements per unit length along each boundary, corresponding to between 25k to 80k total elements depending on domain size as well as number and shapes of voids or inclusions. We employ bilinear quadrilateral CPE4 plain-strain elements for the cases involving a linear elastic material, and their hybrid constant-pressure counterpart CPE4H for the cases involving a hyperelastic material. We apply a load \(P_{o}/E=0.01\) for the cases involving a linear elastic material, and a load \(P_{o}/E=0.173\) for the cases involving a hyperelastic material.

## References

* [1] D. Adalsteinsson and J. A. Sethian. A fast level set method for propagating interfaces. _Journal of computational physics_, 118(2):269-277, 1995.
* [2] A. Adler and D. Holder. _Electrical impedance tomography: methods, history and applications_. CRC Press, 2021.
* [3] G. Allaire and R. Kohn. Topology optimization and optimal shape design using homogenization. In _Topology design of structures_, pages 207-218. Springer, 1993.
* [4] G. Allaire, F. Jouve, and A.-M. Toader. Structural optimization using sensitivity analysis and a level-set method. _Journal of computational physics_, 194(1):363-393, 2004.
* [5] H. B. Ameur, M. Burger, and B. Hackl. Level set methods for geometric inverse problems in linear elasticity. _Inverse Problems_, 20(3):673, 2004.
* [6] H. T. Banks, F. Kojima, and W. P. Winfree. Boundary estimation problems arising in thermal tomography. _Inverse problems_, 6(6):897, 1990.
* [7] M. P. Bendsoe. Optimal shape design as a material distribution problem. _Structural optimization_, 1(4):193-202, 1989.
* [8] M. P. Bendsoe and O. Sigmund. _Topology optimization: theory, methods, and applications_. Springer Science & Business Media, 2003.
* [9] M. Bonnet and A. Constantinescu. Inverse problems in elasticity. _Inverse problems_, 21(2):R1, 2005.
* [10] T. F. Chan and X.-C. Tai. Level set and total variation regularization for elliptic inverse problems with discontinuous coefficients. _Journal of Computational Physics_, 193(1):40-66, 2004.
* [11] Y. Chen and L. Dal Negro. Physics-informed neural networks for imaging and parameter retrieval of photonic nanostructures from near-field data. _APL Photonics_, 7(1):010802, 2022.
* [12] M. Cheney, D. Isaacson, and J. C. Newell. Electrical impedance tomography. _SIAM review_, 41(1):85-101, 1999.
* [13] V. Cherepenin, A. Karpov, A. Korjenevsky, V. Kornienko, A. Mazaletskaya, D. Mazourov, and D. Meister. A 3d electrical impedance tomography (eit) system for breast cancer detection. _Physiological measurement_, 22(1):9, 2001.
* [14] D. Colton, J. Coyle, and P. Monk. Recent developments in inverse acoustic scattering theory. _Siam Review_, 42(3):369-414, 2000.
* [15] A. Daw, J. Bu, S. Wang, P. Perdikaris, and A. Karpatne. Rethinking the importance of sampling in physics-informed neural networks. _arXiv preprint arXiv:2207.02338_, 2022.
* [16] M. Dissanayake and N. Phan-Thien. Neural-network-based approximations for solving partial differential equations. _Communications in Numerical Methods in Engineering_, 10(3):195-201, 1994.
* [17] S. Dong and N. Ni. A method for representing periodic functions and enforcing exactly periodic boundary conditions with deep neural networks. _Journal of Computational Physics_, 435:110242, 2021.
* [18] O. Dorn and D. Lesselier. Level set methods for inverse scattering. _Inverse Problems_, 22(4):R67, 2006.
* [19] O. Dorn, E. L. Miller, and C. M. Rappaport. A shape reconstruction method for electromagnetic tomography using adjoint fields and level sets. _Inverse problems_, 16(5):1119, 2000.
* [20] H. A. Eschenauer, V. V. Kobelev, and A. Schumacher. Bubble method for topology and shape optimization of structures. _Structural optimization_, 8(1):42-51, 1994.
* [21] R. J. Gladstone, M. A. Nabian, and H. Meidani. Fopinns: A first-order formulation for physics informed neural networks. _arXiv preprint arXiv:2210.14320_, 2022.
* [22] H. Griffiths. Magnetic induction tomography. _Measurement science and technology_, 12(8):1126, 2001.
* [23] A. Gropp, L. Yariv, N. Haim, M. Atzmon, and Y. Lipman. Implicit geometric regularization for learning shapes. In _International Conference on Machine Learning_, pages 3789-3799. PMLR, 2020.
* [24] X. Guo, W. Zhang, and W. Zhong. Doing topology optimization explicitly and geometrically?a new moving morphable components based framework. _Journal of Applied Mechanics_, 81(8), 2014.
* [25] E. Haghighat, M. Raissi, A. Moure, H. Gomez, and R. Juanes. A physics-informed deep learning framework for inversion and surrogate modeling in solid mechanics. _Computer Methods in Applied Mechanics and Engineering_, 379:113741, 2021.
* [26] A. Harandi, A. Moeineddin, M. Kaliske, S. Reese, and S. Rezaei. Mixed formulation of physics-informed neural networks for thermo-mechanically coupled systems and heterogeneous domains. _arXiv preprint arXiv:2302.04954_, 2023.
* [27] J. Haslinger and R. A. Makinen. _Introduction to shape optimization: theory, approximation, and computation_. SIAM, 2003.
* [28] C. J. Hellier. _Handbook of nondestructive evaluation_. McGraw-Hill Education, 2013.
* [29] A. Henkes, H. Wessels, and R. Mahnken. Physics informed neural networks for continuum micromechanics. _Computer Methods in Applied Mechanics and Engineering_, 393:114790, 2022.
* [30] M. Kadic, G. W. Milton, M. van Hecke, and M. Wegener. 3d metamaterials. _Nature Reviews Physics_, 1(3):198-210, 2019.
* [31] D. P. Kingma and J. Ba. Adam: A method for stochastic optimization. _arXiv preprint arXiv:1412.6980_, 2014.
* [32] H. T. Kollmann, D. W. Abueidda, S. Koric, E. Guleryuz, and N. A. Sobh. Deep learning for topology optimization of 2d metamaterials. _Materials & Design_, 196:109098, 2020.
* [33] I. E. Lagaris, A. Likas, and D. I. Fotiadis. Artificial neural networks for solving ordinary and partial differential equations. _IEEE transactions on neural networks_, 9(5):987-1000, 1998.
* [34] H. S. Lee, C. J. Park, and H. W. Park. Identification of geometric shapes and material properties of inclusions in two-dimensional finite bodies by boundary parameterization. _Computer Methods in Applied Mechanics and Engineering_, 181(1-3):1-20, 2000.
* [35] L. Lu, X. Meng, Z. Mao, and G. E. Karniadakis. Deepxde: A deep learning library for solving differential equations. _SIAM Review_, 63(1):208-228, 2021.
* [36] L. Lu, R. Pestourie, W. Yao, Z. Wang, F. Verdugo, and S. G. Johnson. Physics-informed neural networks with hard constraints for inverse design. _SIAM Journal on Scientific Computing_, 43(6):B1105-B1132, 2021.
* [37] L. Ma and M. Soleimani. Magnetic induction tomography methods and applications: A review. _Measurement Science and Technology_, 28(7):072001, 2017.
* [38] W. Ma, Z. Liu, Z. A. Kudyshev, A. Boltasseva, W. Cai, and Y. Liu. Deep learning for the design of photonic structures. _Nature Photonics_, 15(2):77-90, 2021.
* [39] Y. Mei, R. Fulmer, V. Raja, S. Wang, and S. Goenezen. Estimating the non-homogeneous elastic modulus distribution from surface deformations. _International Journal of Solids and Structures_, 83:73-80, 2016.
* [40] Y. Mei, Z. Du, D. Zhao, W. Zhang, C. Liu, and X. Guo. Moving morphable inclusion approach: an explicit framework to solve inverse problem in elasticity. _Journal of Applied Mechanics_, 88(4), 2021.
* [41] S. Molesky, Z. Lin, A. Y. Piggott, W. Jin, J. Vuckovic, and A. W. Rodriguez. Inverse design in nanophotonics. _Nature Photonics_, 12(11):659-670, 2018.
* [42] S. Mowlavi and S. Nabi. Optimal control of PDEs using physics-informed neural networks. _Journal of Computational Physics_, 473:111731, 2023.
* [43] S. Osher and R. Fedkiw. _Level set methods and dynamic implicit surfaces_. Springer-Verlag, 2003.
* [44] S. Osher and J. A. Sethian. Fronts propagating with curvature-dependent speed: Algorithms based on hamilton-jacobi formulations. _Journal of computational physics_, 79(1):12-49, 1988.
* [45] M. Raissi, P. Perdikaris, and G. E. Karniadakis. Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. _Journal of Computational physics_, 378:686-707, 2019.
* [46] M. Raissi, A. Yazdani, and G. E. Karniadakis. Hidden fluid mechanics: Learning velocity and pressure fields from flow visualizations. _Science_, 367(6481):1026-1030, 2020.
* [47] C. Rao, H. Sun, and Y. Liu. Physics-informed deep learning for computational elastodynamics without labeled data. _Journal of Engineering Mechanics_, 147(8):04021043, 2021.
* [48] S. Rezaei, A. Harandi, A. Moeineddin, B.-X. Xu, and S. Reese. A mixed formulation for physics-informed neural networks as a potential solver for engineering problems in heterogeneous domains: comparison with finite element method. _Computer Methods in Applied Mechanics and Engineering_, 401:115616, 2022.
* [49] L. Robledo, M. Carrasco, and D. Mery. A survey of land mine detection technology. _International Journal of Remote Sensing_, 30(9):2399-2410, 2009.
* [50] F. Sahli Costabal, Y. Yang, P. Perdikaris, D. E. Hurtado, and E. Kuhl. Physics-informed neural networks for cardiac activation mapping. _Frontiers in Physics_, 8:42, 2020.
* [51] O. Sigmund and K. Maute. Topology optimization approaches. _Structural and Multidisciplinary Optimization_, 48(6):1031-1055, 2013.
* [52] V. Sitzmann, J. Martel, A. Bergman, D. Lindell, and G. Wetzstein. Implicit neural representations with periodic activation functions. _Advances in Neural Information Processing Systems_, 33:7462-7473, 2020.
* [53] N. Sukumar and A. Srivastava. Exact imposition of boundary conditions with distance functions in physics-informed deep neural networks. _Computer Methods in Applied Mechanics and Engineering_, 389:114333, 2022.
* [54] C. L. Wight and J. Zhao. Solving Allen-Cahn and Cahn-Hilliard equations using the adaptive physics informed neural networks. _Communications in Computational Physics_, 29(3):930-954, 2021.
* [55] J. C. Wong, C. Ooi, A. Gupta, and Y.-S. Ong. Learning in sinusoidal spaces with physics-informed neural networks. _arXiv preprint arXiv:2109.09338_, 2021.
* [56] E. Zhang, M. Dao, G. E. Karniadakis, and S. Suresh. Analyses of internal structures and defects in materials using physics-informed neural networks. _Science advances_, 8(7):eabk0644, 2022.


**Figure 3**: **Identification of inclusions in a linear elastic matrix.** Final Cauchy stress components \(\sigma_{xx}\) (**a**), \(\sigma_{yy}\) (**b**), and \(\sigma_{xy}\) (**c**), displayed in the deformed configuration obtained from the final displacement components \(u_{1}\) and \(u_{2}\), for the cases reported in Fig. 4e. The grey dotted lines show the outline of the matrix surface in the reference configuration.
Extended Data Fig. 5. **Effect of sparse measurements on the identification of voids in a linear elastic matrix with one circle-shaped void.** Final material density obtained in our framework when using fewer measurement locations and restricting the number of measurements to a subset of the outer surfaces.

Extended Data Fig. 6. **Effect of sparse measurements on the identification of voids in a linear elastic matrix with one star-shaped and one rectangle-shaped void.** Final material density obtained in our framework when using fewer measurement locations and restricting the number of measurements to a subset of the outer surfaces.
Extended Data Fig. 7. **Effect of sparse measurements on the identification of voids in a linear elastic matrix with one slit-shaped void.** Final material density obtained in our framework when using fewer measurement locations and restricting the number of measurements to a subset of the outer surfaces.
**Supplementary information for:**

**Topology optimization with physics-informed neural networks:**

**application to noninvasive detection of hidden geometries**

Saviz Mowlavi\({}^{1,\,2,}\)1, and Ken Kamrin\({}^{1,\,\dagger}\)2

Footnote 1: mowlavi@merl.com

Footnote 2: kkamrin@mit.edu

\({}^{1}\)_Department of Mechanical Engineering, MIT, Cambridge, MA 02139, USA_

\({}^{2}\)_Mitsubishi Electric Research Laboratories, Cambridge, MA 02139, USA_

## I Uniqueness of solutions

We sketch a proof for the uniqueness of solutions to the geometry detection elasticity problem considered in this paper. For the specific case of a two-dimensional linear elastic material with a single void, it has been proved that there exists at most one cavity which yields the same surface displacements and stresses on an finite portion of the external boundary [1]. Our approach is different and applicable to any three-dimensional problem governed by elliptic PDEs (be they linear or nonlinear, heat conduction, elasticity, etc) and containing multiple voids or rigid inclusions. Before sketching the proof, we state a few lemmas that will be useful:

* **Lemma 1**. Any \(C^{2}\) displacement field that solves the equilibrium Navier-Cauchy equations of linear elasticity with constant elastic moduli is a real analytic function of space. Ths holds because the Navier-Cauchy equations are elliptic, so any \(C^{2}\) solution must be real analytic [3].
* **Lemma 2**. Any real analytic function on a connected domain \(\Omega\) that vanishes on a finite and connected subset \(\mathcal{S}\subset\Omega\) necessarily vanishes in all of \(\Omega\). This is a specialization of the identity theorem for analytic functions [2].
* **Lemma 3**. If \(\mathbf{u}\) is a \(C^{2}\) displacement solution to the Navier-Cauchy equations on a connected domain \(\Omega\), and both \(\mathbf{u}\) and traction \(\mathbf{t}\) vanish on a smooth finite portion of the boundary \(S\subset\partial\Omega\), then \(\mathbf{u}\) must necessarily vanish in all of \(\Omega\). This fact can be seen by extending the domain \(\Omega\) along \(S\) by some amount \(\Omega^{ext}\) and defining \(\mathbf{u}=\mathbf{0}\) there. Then \(\mathbf{u}\) satisfies the Navier-Cauchy equations at each point in \(\Omega\cup\Omega^{ext}\), and it can be shown that all second derivatives exist and are continuous on \(S\). Thus, on \(\Omega\cup\Omega^{ext}\), \(\mathbf{u}\) is an analytic function (by Lemma 1) that vanishes in \(\Omega^{ext}\), and therefore \(\mathbf{u}\) vanishes in all of \(\Omega\cup\Omega^{ext}\) (by Lemma 2) and thus in \(\Omega\).

We sketch the proof for the case of a single void or inclusion, but the same reasoning generalizes to any number of voids or inclusions. Consider two bodies \(\mathcal{B}^{(1)}\) and \(\mathcal{B}^{(2)}\) with the same material properties and sharing the same external boundary \(\partial\mathcal{B}^{\mathrm{ext}}\) (Supplementary Fig. 1a,b). Each body contains a single smooth void or rigid inclusion, characterized by the connected domains \(\mathcal{I}^{(1)}\) and \(\mathcal{I}^{(2)}\) with boundaries \(\partial\mathcal{I}^{(1)}\) and \(\partial\mathcal{I}^{(2)}\), respectively. The two bodies are subjected to the same external loading, which consists of an applied displacement \(\mathbf{u}\) on a portion \(\partial\mathcal{B}^{\mathrm{ext}}_{u}\) of \(\partial\mathcal{B}^{\mathrm{ext}}\) and an applied traction \(\mathbf{\hat{t}}\) on a portion \(\partial\mathcal{B}^{\mathrm{ext}}_{t}\) of \(\partial\mathcal{B}^{\mathrm{ext}}\), with \(\partial\mathcal{B}^{\mathrm{ext}}=\partial\mathcal{B}^{\mathrm{ext}}_{u} \cup\partial\mathcal{B}^{\mathrm{ext}}_{t}\). This loading generates a displacement and stress solution within each body, which we denote by \(\mathbf{u}^{(1)}\), \(\boldsymbol{\sigma}^{(1)}\) and \(\mathbf{u}^{(2)}\), \(\boldsymbol{\sigma}^{(2)}\), respectively. These solutions must be nontrivial, meaning that the surface tractions \(\mathbf{t}^{(1)}\) and \(\mathbf{t}^{(2)}\) do not vanish everywhere on \(\partial\mathcal{B}^{\mathrm{ext}}\). Finally, we assume that we measure identical surface displacements, i.e. \(\mathbf{u}^{(1)}=\mathbf{u}^{(2)}=\mathbf{u}^{m}\), on a finite portion \(\partial\mathcal{B}^{m}\) of \(\partial\mathcal{B}^{\mathrm{ext}}_{t}\).

We will prove that there cannot be two distinct shapes \(\mathcal{I}^{(1)}\) and \(\mathcal{I}^{(2)}\) yielding nontrivial solutions that are identical on \(\partial\mathcal{B}^{m}\). Let us first subtract the displacement and stress solutions of the two bodies, which yields a displacement field \(\Delta\mathbf{u}=\mathbf{u}^{(1)}-\mathbf{u}^{(2)}\) and a stress field \(\Delta\boldsymbol{\sigma}=\boldsymbol{\sigma}^{(1)}-\boldsymbol{\sigma}^{(2)}\) defined over the intersection \(\mathcal{B}^{(1)}\cap\mathcal{B}^{(2)}\) between the two bodies (Supplementary Fig. 1c). This 'difference solution' itself satisfies the governing equations of linear elasticity owing to their linearity, and its displacement and traction vanish on \(\partial\mathcal{B}^{m}\). Therefore, we can show using Lemmas 1 and 3 that the difference solution vanishes in its domain of definition, meaning that \(\mathbf{u}^{(1)}=\mathbf{u}^{(2)}\) and \(\boldsymbol{\sigma}^{(1)}=\boldsymbol{\sigma}^{(2)}\) in \(\mathcal{B}^{(1)}\cap\mathcal{B}^{(2)}\). We now focus on the behavior of the solution in body \(\mathcal{B}^{(1)}\) and treat separately the cases of a void or rigid inclusion.

In the case of a void (Supplementary Fig. 1d), the traction \(\mathbf{t}^{(1)}\) on the internal boundary \(\partial\mathcal{I}^{(1)}\) vanishes. Using the fact that \(\boldsymbol{\sigma}^{(1)}=\boldsymbol{\sigma}^{(2)}\) in \(\mathcal{B}^{(1)}\cap\mathcal{B}^{(2)}\), we know that \(\mathbf{t}^{(1)}\) also vanishes on \(\mathcal{B}^{(1)}\cap\partial\mathcal{I}^{(2)}\). As a result, there exists a finite region \(\tilde{\mathcal{B}}^{(1)}=\mathcal{B}^{(1)}\cap\mathcal{I}^{(2)}\) with zero boundary traction. Thus, \(\boldsymbol{\sigma}^{(1)}\) must vanish not only in \(\tilde{\mathcal{B}}^{(1)}\), but also in the
entire body \(\mathcal{B}^{(1)}\) as a consequence of Lemma 2. This violates the fact that the surface traction \(\mathbf{t}^{(1)}\) cannot vanish everywhere on \(\partial\mathcal{B}^{\mathrm{ext}}\).

In the case of a rigid inclusion (Supplementary Fig. 1e), the displacement \(\mathbf{u}^{(1)}\) on the internal boundary \(\partial\mathcal{I}^{(1)}\) corresponds to that of a rigid motion, i.e. \(\mathbf{u}^{(1)}=\mathbf{u}^{(1)}_{r}+\mathbf{\theta}^{(1)}_{r}\times\mathbf{x}\) for some fixed \(\mathbf{u}^{(1)}_{r}\) and \(\mathbf{\theta}^{(1)}_{r}\). Using the fact that \(\mathbf{u}^{(1)}=\mathbf{u}^{(2)}\) in \(\mathcal{B}^{(1)}\cap\mathcal{B}^{(2)}\), we know that \(\mathbf{u}^{(1)}=\mathbf{u}^{(2)}=\mathbf{u}^{(2)}_{r}+\mathbf{\theta}^{(2)}_{r} \times\mathbf{x}\) on \(\mathcal{B}^{(1)}\cap\partial\mathcal{I}^{(2)}\) for some other fixed \(\mathbf{u}^{(2)}_{r}\) and \(\mathbf{\theta}^{(2)}_{r}\). However, given that these two rigid motions must coincide at \(\partial\mathcal{I}^{(1)}\cap\partial\mathcal{I}^{(2)}\), we must have \(\mathbf{u}^{(1)}_{r}=\mathbf{u}^{(2)}_{r}\) and \(\mathbf{\theta}^{(1)}_{r}=\mathbf{\theta}^{(2)}_{r}\). As a result, the entire surface of the finite region \(\tilde{\mathcal{B}}^{(1)}=\mathcal{B}^{(1)}\cap\mathcal{I}^{(2)}\) undergoes a rigid motion. Thus, \(\mathbf{\sigma}^{(1)}\) must vanish not only in \(\tilde{\mathcal{B}}^{(1)}\), but also in the entire body \(\mathcal{B}^{(1)}\) as a consequence of Lemma 2. As in the case of the void, this violates the fact that the surface traction \(\mathbf{t}^{(1)}\) cannot vanish everywhere on \(\partial\mathcal{B}^{\mathrm{ext}}\).
## IIIL-Posed Nonlinear Thermal Problem

To illustrate the flexibility of our TO framework with respect to the physical model and the sparsity of the data, in this appendix we apply our method to an ill-posed nonlinear thermal imaging problem. Consider a nonlinearly conducting matrix (Supplementary Fig. 2a) whose heat flux \(\mathbf{q}\) relates to the temperature \(T\) through a nonlinear Fourier's law \(\mathbf{q}=-k(1+T/T_{0})\nabla T\), where \(k\) is akin to a thermal conductivity and \(T_{0}\) is a reference temperature. The matrix contains a hidden inclusion that is either perfectly insulating (no inside heat flux) or perfectly conducting (no inside temperature gradient). A thermal loading is applied, consisting of a prescribed unit temperature on the left boundary, zero temperature on the right boundary, and insulated top and bottom boundaries (Supplementary Fig. 2a). The goal of the inverse problem is to identify the location and shape of the inclusion under different ill-posed scenarios where the applied boundary condition and the measurements are both assumed unavailable on one of the four sides, while the temperature or flux resulting from the prescribed thermal loading are measured on the remaining three sides (Supplementary Fig. 2b).

The experiment is simulated in the FEM software Abaqus using biquadratic DC2D8 diffusive heat transfer elements, and considering \(k=1\) and \(T_{0}=1\). The inverse problem is solved with our PINN-based TO framework by constructing neural network approximations for the physical quantities \(\mathbf{\psi}=(T,\mathbf{q})\) and the density field \(\rho\). Similarly to the elasticity examples, the neural networks are designed to inherently satisfy the boundary conditions (on the three sides where they are assumed to be known). The governing equations included in the loss function comprise the conservation law \(\nabla\cdot\mathbf{q}=0\) as well as the nonlinear Fourier's law \(F(\mathbf{q},T,\rho)=0\). The latter, expressed over both the matrix and the inclusion, takes the form \(\mathbf{q}+\rho k(1+T/T_{0})\nabla T=0\) in the presence of a perfectly-insulating inclusion, or \(\rho\mathbf{q}/k+(1+T/T_{0})\nabla T=0\) in the presence of a perfectly-conducting inclusion. We use the same training parameters and the same weights in the loss function as in the square elastic matrix examples.

Despite the complete lack of information along an entire side, with the applied boundary condition and measurements both missing, our TO framework is able to detect the slit-shaped inclusion with reasonable accuracy in both the perfectly-insulating case (Supplementary Fig. 2c) and the perfectly-conducting case (Supplementary Fig. 2d). This example showcases the ability of the framework to generate good results when even the forward problem is ill-posed.
