{'fraction_non_alphanumeric': 0.03078543484802889, 'fraction_numerical': 0.014234125789948841, 'mean_word_length': 4.8279551034724655, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
SELTO: Sample-Efficient Learned Topology Optimization

We present a sample-efficient deep learning strategy for topology optimization. Our end-to-end approach is supervised and includes physics-based preprocessing and equivariant networks. We analyze how different components of our deep learning pipeline influence the number of required training samples via a large-scale comparison. The results demonstrate that including physical concepts not only drastically improves the sample efficiency but also the predictions' physical correctness. Finally, we publish two topology optimization datasets containing problems and corresponding ground truth solutions. We are confident that these datasets will improve comparability and future progress in the field.

1. Introduction. The field of topology optimization (TO) is a computational discipline aiming to optimize mechanical structures that fulfill user-specified constraints and minimize user-defined cost functions. Since its development in 1988 [4], TO has become a powerful tool adopted by engineers in many fields, ranging from heat transfer [15] and acoustics [17,49] to fluid [8] and solid mechanics [19].
Solid Isotropic Material with Penalization (SIMP) [5], arguably the most important classical TO method, uses a density-based setup. Here one discretizes a given design domain into voxels, each having a density value between 0 and 1. These values represent the degree of material present in the voxel and thereby define a mechanical structure. One employs the partial differential equation (PDE) governing the physical phenomenon of interest -we will focus on linear elasticity -to determine the performance regarding the specified constraints and cost function. One then uses iterative optimization methods to adjust the voxel densities to improve the performance. From a mathematical perspective, this constitutes a PDE parameter identification problem.
Density-based methods work well for low-resolution voxel meshes. However, due to their iterative nature, methods like SIMP repeatedly solve the PDE of the physical phenomenon of interest. Therefore, increasing the resolution of the voxel mesh becomes progressively computationally demanding to the point of practical impossibility [1]. Recent research is trying to speed up and improve this optimization process using modern machine learning algorithms, most notably deep learning (DL), i.e., neural networks [6]. Still, the subdiscipline of applying DL to TO is in its infancy, with most papers focusing on two-dimensional settings [27,31,32,36,42,45,52,53].
Applying DL to TO has its advantages and disadvantages. Unlike classical methods, many DL models can solve TO problems on a subsecond-scale. Critically though, this assumes the existence of a trained model. This constitutes a considerable impediment as many DL models only function on load and boundary conditions similar to their training data -the so-called generalization problem [25]. Generalization would be less critical if training a DL model would not require large-scale training datasets, providing TO problems and their respective optimal structures.
We want to address these issues via the following contributions: • We show that enforcing rotational and mirror equivariance in DL models improves generalization and significantly reduces the dependence on large training sets. We believe this novel incorporation of fundamental physical laws to be a significant step toward real-world applicability. • We study sample efficiency via a large-scale comparison of different DL architectures and input preprocessing strategies. • We provide two datasets containing problems and corresponding ground truth solutions, hoping to improve the comparability in the field.
2. Preliminaries and related work. We now briefly introduce classical TO, followed by a review of the current literature on DL for TO.
2.1. Density-based topology optimization. We start by discussing densitybased TO, arguably the most common classical TO approach -also used to generate our datasets (see Section 3).
Density-based TO [5] aims to minimize a cost or objective function by adjusting the material's density distribution ρ : Ω → {ρ min , 1} over a fixed domain Ω ⊂ R d , typically d ∈ {2, 3}. Here 0 < ρ min 1 is the minimal density, where ρ min > 0 is a numerical necessity for solving the governing PDE. Additionally, this optimization is subject to physical constraints. The specified objective function and constraints may vary depending on the user's needs.
This paper focuses on compliance minimization, the most common setting for mechanical problems. The corresponding optimization problem reads as follows: Here (2.1a) is the compliance objective function, F represents the global load distribution, u are the displacement and K is the symmetric positive operator of linear elasticity. K includes the characteristic properties of the used material described by Young's modulus and Poisson's ratio, which we denote by E ∈ R and ν ∈ [0, 0.5], respectively. One constrains the amount of allowed material by V max (2.1c), and one may include a stress constraint (2.1d) to ensure that the maximal von Mises stress σ vM is below the yield stress σ yield of the material.
In practice, the SIMP [5] method relaxes the problem. Instead of strictly binary density values, one allows ρ : Ω → [ρ min , 1] and encourages near-binary densities by extending Young's modulus over ρ's interval via E(ρ) = E 0 ρ p , where E 0 is the original (isotropic) material's modulus. The exponent p controls the penalization of non-binary densities, usually p = 3. For implementation purposes, one discretizes the design space Ω into a regular voxel grid and iteratively updates ρ until a user-specified convergence criterion is met.
As discussed in the introduction, despite several advancements in structural TO, one of the main challenges is the high computational cost. The main bottleneck in classical iterative approaches is that each iteration uses the displacements and stresses for the current density and therefore has to solve the PDE for linear elasticity in (2.1b). Due to this computational challenge, TO at high resolutions, i.e., over many voxels, can take hours, even days [1]. This has inspired researchers to develop DL-based TO methods to reduce or even eliminate the need to solve PDEs.

2.2.
Neural networks for topology optimization. One can broadly classify advances in TO using DL into three categories [53]: Reduce or eliminate SIMP iterations: The first serious attempts of TO via DL aimed to reduce the number of SIMP iterations [3,42,48]. In 2017, [42] interpreted two-dimensional TO problems as image-to-image regression problems and were the first to apply convolutional neural networks (CNNs) to TO. Following wellknown image processing approaches, they trained a UNet model [40] to map from intermediate SIMP iterations to the final structure. In 2018, [3] transferred these ideas to the three-dimensional case. The authors of [48] made each SIMP iteration cheaper by running it on a coarse resolution. They then applied a DL-based superresolution method to increase the structure's final granularity. In 2020, [2] were the first to use residual neural networks (ResNets) [54] and to consider two-dimensional nonlinear elasticity.
In 2019, [50] developed the first end-to-end learning routine that directly predicts the final density without performing any SIMP iterations. They also create a generative framework to increase the resolution of their predicted designs. Based on that, [31,32,39,41] conducted a series of explorations on generative adversarial network (GAN)-based TO algorithms [21].
Compared to previous research, [31] and [52] achieved a better generalization by not directly giving the network the boundary conditions as input. Instead, they passed displacements and von Mises stresses into the network. They argue that neural networks have difficulties extending to previously unseen boundary conditions if the input data is very sparse since the high sparsity of the input matrices leads to high variance of the mapping function.
Substitute SIMP's PDE solver: These methods aim to remove classical PDE solvers from the SIMP algorithm, removing its primary bottleneck. [36] proposes a dual-model neural network using a forward model to compute the compliance of the structure and an adjoint model to determine the derivatives with respect to the density of each voxel. Similarly, [12] and [27] use neural networks to replace the gradient and objective function computation.
Neural reparameterization: Using implicit neural representation for complex signals is an ongoing research topic in computer vision [11] and engineering [29]. Several TO papers [16,23,51,53] use neural networks to reparameterize the TO's density field. While [23] maps latent vectors to discrete grid densities, [10] uses multilayer perceptrons to learn a continuous mapping from spatial locations to density values. Since these models are mesh-independent, they can represent the density function at arbitrary resolutions. However, one usually still requires PDE evaluations for training, which is computationally demanding.
While most DL approaches to TO are entirely oblivious to the underlying physics, there is a small number of approaches that have started to include physics-inspired properties into the training process: In [3,37], the authors augment the dataset by including rotations and mirrors of given loads and boundary conditions to encourage equivariance. [31] and [54] involve physics by feeding the network strain and stress information as inputs. In [9] and [53], the authors introduce physics via the loss function design, comparable to physics-informed neural networks (PINNs) [38].
To our knowledge, no literature incorporated the problem's underlying physics by modifying the architecture of the DL model itself. We demonstrate that we can dramatically improve sample efficiency, i.e., the model's performance when trained on a few training samples, and facilitate geometric reasoning by restricting the hypothesis space to group equivariant models. [14] introduced the first group equivariant CNNs in 2016. Nowadays, their application ranges from chemistry [46] and physics [7] to a wide range of tasks in geometric DL [44,20]. For image processing tasks, [18] showed that enforcing relevant equivariances can improve generalization performance and robustness.

Datasets.
We find a substantial lack in the availability of public threedimensional TO datasets, i.e., to the best of our knowledge there is none. This is problematic for several reasons. First, it is hard to reproduce other published results. Second, it requires each researcher to generate their own datasets, which is time-consuming and computationally demanding. Lastly, the lack of established TO datasets impedes the comparability of results throughout the community. To alleviate this issue, we publish two three-dimensional TO datasets containing samples of mounting brackets, which we call disc dataset and sphere dataset 1 .
Each dataset consists of TO problems and associated ground truth density distributions. We generated the samples in cooperation with the ArianeGroup 2 and ELISE 3 using the Altair OptiStruct 4 implementation of SIMP within the ELISE software. The ArianeGroup designed the mounting brackets in the datasets to be of practical use, though real-world aerospace applications would require more complex load cases. The samples are discretized on a n x × n y × n z voxel grid, where the choice of n x , n y , n z ∈ N varies depending on the dataset. Both datasets have fixed Dirichlet boundary conditions but variable force positions and magnitudes.
For convenience we introduce the following definition: Definition 3.1. Let k ∈ N and X ⊆ R. If x ∈ X k×nx×ny×nz , then we call x a k-channel X-tensor.
One can uniquely characterize each TO problem via the following properties: 1. The number of voxels (n x , n y , n z ) and the voxel size in millimeters in each direction. 2. Material properties, given by Young's modulus E, Poisson's ratio ν and a yield stress criterion σ yield . We choose E = 70 GPa, ν = 0.3 and σ yield = 450 MPa for both datasets. 3. A 3-channel {0, 1}-tensor ω Dirichlet to encode the presence of directional Dirichlet boundary conditions for every voxel. 1s indicate the presence, and 0s the absence of homogeneous Dirichlet boundary conditions. 4. A 3-channel R-tensor F to encode external forces, given in N/m 3 . The three channels describe the force magnitudes in each spacial dimension. 5. A 1-channel {−1, 0, 1}-tensor ω design to encode design space information. We use 0s and 1s to constrain voxel densities to be 0 or 1, respectively. Entries of −1s indicate a lack of density constraints, which signifies that the density in that voxel can be freely optimized. This naturally defines the voxel sets Ω −1 , Ω 0 , and Ω 1 . For voxels that have Dirichlet boundary conditions or loads assigned to them we enforce the density value to be 1 by setting ω design = 1.
Note that all tensors are defined voxel-wise, including ω Dirichlet and F . This should make our datasets easy to use in DL applications since it allows for a shape-consistent tensor representation. The SIMP method does not always provide a physically plausible solution for a TO problem, i.e., some solutions break under their load cases. Therefore, we clean both datasets after the dataset generation process by rejecting failing samples. In our case, this decreased dataset sizes by around 20 -25%.

Disc dataset.
The disc dataset consists of 9328 samples, each having 39 × 39 × 4 voxels. The dataset comprises 1709 samples with a load case with a single point of attack -we call those the disc simple dataset -and 7619 samples with load cases having two points of attack each -we call those the disc complex dataset. We split the disc simple and disc complex dataset into training and validation subsets, with 200 validation samples each. We refer to the combination of the disc simple and disc complex dataset as the disc combined dataset. We present ground truth examples from the disc dataset in Figures 1a and 1b. Due to the reduced number of voxels in the z-direction, this dataset is relatively easy to handle while not being a two-dimensional dataset.
3.2. Sphere dataset. The sphere dataset consists of 604 samples, 188 in sphere simple, and 416 in sphere complex. Each sample has 39 × 39 × 21 voxels. We split each dataset into training and validation subsets with 38 validation samples. This ensures at least 150 training samples per dataset. Mutatis mutandis to disc we define the sphere combined dataset. We present exemplary ground truths in Figures 1c and 1d.

Methods.
This section introduces the different DL pipelines we will evaluate and analyze in Section 5. We examine different input preprocessing strategies and the effects of physically motivated group averaging [34].
We use an end-to-end learning approach, i.e., train a neural network f to map preprocessed TO problems to an optimized density distribution. Since both the disc and the sphere dataset contain ground truth densities for each problem, we can train f in a supervised fashion.

Preprocessings.
Choosing a suitable input preprocessing strategy is crucial for DL. We now present the two main preprocessings we use in this paper. It is possible to combine these via simple channel-wise concatenation of their outputs.
1. Trivial preprocessing. The input of the neural network is a 7-channel tensor which results from the channel-wise concatenation of Dirichlet boundary conditions ω Dirichlet , loads F , and design space information ω design . Additionally, we normalize each sample's F via the mean F ∞ over all training samples. This is arguably the most straightforward and intuitive type of preprocessing for learned end-to-end TO. 2. PDE preprocessing. As proposed in [52], we first define an initial density distribution ρ init that is 1 on Ω 1 and Ω −1 , and 0 on Ω 0 . For ρ init , we then compute the initial von Mises stresses, which we obtain by solving the PDE for linear elasticity. We normalize the resulting tensor analogously to the normalization of F above. These initial von Mises stresses are then used as a 1-channel input to the neural network. Analogously, it would also be possible to use the full initial stress tensor or the initial displacements as network input.
We illustrate both preprocessing strategies and our model pipeline in Figure 2.

Architecture.
We choose a UNet [40] as our neural network architecture, which is a convolutional encoder-decoder network. The encoder consists of repeated application of convolutions, each followed by a rectified linear unit (ReLU) activation function and a max pooling operation. During the encoding, the encoder of the UNet reduces the spatial information while it increases the feature information. The decoder then extends the feature and spatial information through convolution and upsampling steps and concatenations with high-resolution features from the encoder. See Appendix A for more details on the UNet architecture. Fig. 2: Illustration of our model pipeline, using an end-to-end learning approach, where the function f represents the neural network. Top: Trivial preprocessing. We concatenate ω Dirichlet , ω design , and forces F = F 1 ∪ F 2 into a 7-channel tensor, which we use as input for f . Bottom: PDE preprocessing. We compute an initial density distribution ρ init and corresponding initial displacements and stresses. The image displays the initial von Mises stress distribution, which is a 1-channel tensor.

Equivariance.
Equivariance is the property of a function to commute with the actions of a symmetry group acting on its domain and range. For a given transformation group G, we say that a function f : where T X g and T Y g denote linear group actions in the corresponding spaces X and Y [14]. That is, transforming an input x by a transformation g and then passing it through f should give the same result as first mapping x through f and then transforming the output f (x) (see left image of Figure 3). In many machine learning tasks, we possess prior knowledge about equivariances our predictor should have. Including such knowledge directly into the model can significantly facilitate learning by freeing up model capacity for other factors of variation [47]. Since G is a group, it also contains the identity transformation and a unique inverse transformation g −1 for each g ∈ G. Therefore, we can reformulate (4.1) as allowing the implementation of equivariance via group averaging [30,34] by defining an equivariance wrapper F f G as Here, g is represented by a 90°clockwise rotation. Right: Illustration of the equivariance wrapper F f G . We transform the input x via group actions T X gi and feed the transformed data into the model f . We then reverse the initial transformation on the outputs via T Y g −1 i and aggregate the results. This process is called group averaging. The choice of the (finite) transformation group G is entirely free.
For an illustration of F f G see the right image of Figure 3. The most popular approach to enforcing group equivariance in DL is enforcing invariance of the network's convolutional filters and activation functions [14]. However, performing group averaging has several advantages. First, it is straightforward to implement. Second, its plug-and-play nature allows the trivial application to any (finite) transformation group G and any model f . In particular, the applicability to any f allows for direct comparisons between any non-equivariant model and its equivariant counterpart.
Natural choices of transformation groups for our applications are the dihedral symmetry group D 4 and the octahedral symmetry group O h , which describe all combinations of 90°rotations and reflections in R 2 and R 3 , respectively. Thus the equivariance wrapper applies 8 transformations in the D 4 and 48 in the O h case. Note that group actions on vector fields do -unlike on scalar fields -affect spacial locations and vector directions. The necessity to account for vector directions makes equivariance wrappers especially suited and flexible to employ.
We should clarify that the equivariance wrapper fundamentally differs from data augmentation, i.e., augmenting the dataset with rotations and reflections. Augmentation averages over the objective, whereas the wrapper averages over the model. This difference endows the wrapper with several significant advantages, most notably: The model does not need to learn about equivariance as it is hardcoded. Consequently, the wrapper is not only approximately equivariant but exactly -importantly, this holds for any input, even for inputs drastically different from the training data.

Numerical experiments.
In this section, we conduct several numerical experiments to illustrate the effectiveness of adding physics-based information via PDE preprocessing and equivariance to our DL models.

Training.
We train and compare different combinations of preprocessings described in Section 4, each with and without equivariance. Due to the reduced number of voxels in the z-direction and for simplicity, we chose the dihedral symmetry group, D 4 , as the transformation group for all our experiments. All models are implemented in PyTorch [33] and trained on an Nvidia GeForce GTX 1080 Ti GPU with 11GB VRAM. We determine the batch sizes individually, depending on the memory capacity and comparability (see Appendix A). We choose the weighted binary cross-entropy (BCE) as our loss function, and calculate the weighting factor based on the training dataset. We use the Adam optimizer [26] with a learning rate of 10 −3 in all our experiments.
We train all models until their improvement on the validation set stalled for 100 epochs and then pick the model corresponding to the best validation epoch; this stopping criterion results in models trained for 30 to 970 epochs.

Evaluation.
We now discuss our evaluation methodology for the analysis of our DL models. In particular, we want to compare the impact of different preprocessing strategies and equivariances on our models. We begin by defining our evaluation criteria.

Evaluation criteria.
We now introduce the criteria we use to evaluate our models over the validation datasets. Before we apply the criteria we first binarize the densities produced by the DL models, such that they only contain 0s and 1s. We make use of the following two criteria: • IoU: In contrast to our loss, the BCE, which is distribution-based, the Intersection over Union (IoU) is region-based. It is defined as where TP, FN and FP denote the number of true positive, false negative and false positive voxel predictions. We limit the evaluation of the IoU to the editable design space Ω −1 . Following Goodhart's law [22] and established practices in the segmentation community, we use the IoU as our primary evaluation metric [28,40]. • Fail percentage: We consider a prediction failed if the von Mises stress in any voxel exceeds the yield stress by more than 10% or if any voxel with a load case does not connect to one containing Dirichlet conditions. The fail percentage is the fraction of failed parts. This criterion makes us the first in the DL for TO community who verify their predictions' mechanical integrity.
A typical evaluation criterion in the DL for TO literature is the (balanced) accuracy metric. However, we found that in most cases, IoU yields comparable values but, in general, reflects the quality of its input more appropriately when used across different datasets. Given that we can interpret the binary densities of our structures as segmentation masks, IoU's applicability is not surprising as it is the most common metric for semantic segmentation. Therefore we will use IoU as our main evaluation criterion.

Sample efficiency.
One of the main impediments to learned supervised TO is the high cost associated with generating new datasets. The necessity to generate thousands of problems with corresponding ground truth solutions often hinders practical applicability in real-world situations since the generation of large datasets can take days, even weeks. Consequently, reducing the number of required training samples is highly beneficial, e.g., by modifying the DL model design. Therefore, we put a particular emphasis on the visualization and measurement of our models' sample efficiency, i.e., the model's performance when trained on few training samples.
We visualize the sample efficiency of a model with what we call sample efficiency curves (SE curves). For each SE curve, we train separate instances of a given model setup on subsets of the original training dataset of varying sizes. We then evaluate and compare the performance of these models on a fixed validation dataset, using IoU and fail percentage as our evaluation criteria. For the disc dataset, we choose training subsets of sizes 2, 4, 10, 50, 100, 150, 250, 500, 1000, and 1500. For the sphere dataset we train on 2, 4, 10, 50, 100, and 150 samples. This way, we obtain an individual SE curve for each evaluation criterion. In order to compare different SE curves of the same criterion, we report two metrics: 1. The normalized area under the curve (AUC) of that criterion up to a training sample size of 150, denoted by AUC 150 . We use this metric to quantify sample efficiency. 2. The final score, which is the value of that criterion achieved by the model trained on the largest training subset.

Results.
This section gives an overview of our numerical results and compares the performance of different models. We begin with our main results, followed by an analysis of our models' generalization capabilities. Finally, we discuss model alterations that did not yield improvements but might still be helpful for further understanding and research.

Main results.
As expected, model performance tends to improve with increased training data. We also observe dramatic boosts in the UNet's performance when incorporating physics via equivariance and trivial+PDE preprocessing. These improvements are especially visible for low numbers of training samples; see the SE curves in Figure 4 and AUC 150 scores in Table 2. From Figure 4, we observe that using trivial+PDE preprocessing and equivariance in combination leads to a reduction in the required training samples by two orders of magnitude while maintaining equivalent IoU scores. Additionally, we notice a reduction in the fail percentage below 20%, comparable to the fail percentage of the classical SIMP algorithm during the dataset creation. This is a remarkable improvement in the predictions' mechanical integrity.
The improvements are also evident when visually comparing ground truths with predictions; see Tables 4 to 7. Further, we observe that adding PDE preprocessing and equivariance leads to predictions closer to the ground truth densities, which is especially noticeable on small training sets. Moreover, applying our equivariance  wrapper does not significantly increase training duration, see Table 8 in Appendix A.
For an unbiased impression of our model predictions and to avoid cherry-picking, we display 20 random samples and predictions from both datasets in Tables 9 and 10 in Appendix C. We show two predictions for each sample, one using trivial preprocessing only and one using trival+PDE preprocessing and equivariance. We trained the UNets on 1500 and 150 samples for disc and sphere, respectively.

Generalization.
As in all machine learning, generalization is also a problem in its application to TO. The problem led to some publications addressing the issue [31,52]. We quantify our models' generalizability by cross-evaluating them on different data subsets not used during training, i.e., we evaluate all models trained on either disc simple, disc complex, or disc combined on each of the others. We proceed analogously for the sphere models and datasets.
We compare the AUC 150 and final scores for each model based on the corresponding SE curve. For both IoU and fail percentage, we present the results for the AUC 150 score in Table 2 and the final score in Table 3. We observe that the addition of PDE preprocessing and equivariance improves the generalization capabilities of our models considerably for both IoU and fail percentage, especially on small training sets. See Table 1, for examples of cross-evaluation predictions on disc and sphere.

Ablations and negative results.
This section discusses several model alterations that did not lead to noticeable improvements in our numerical experiments. We think these results are still valuable to the community.
Above, we exclusively used the dihedral symmetry group, D 4 , in the equivariance wrapper. For the sphere dataset, we also examined the performance of the octahedral symmetry group O h . While still constituting a significant benefit over not utilizing equivariance, this method was inferior to the use of the D 4 symmetry group. We suspect this is due to the fixed location of the sphere dataset's Dirichlet condition, negating the benefit of generalizing to more diverse locations of boundary conditions.
In addition to the preprocessings we examined in Subsection 5.3.1, we experimented with various other preprocessing combinations. We found that convex hull preprocessing (see Appendix B) did not lead to any measurable improvements on the validation dataset, neither on its own nor in combination with other preprocessings. We speculate that this is due to the convex hull being invariant to the direction of the applied forces. We also found that applying PDE preprocessing on its own produced similar but somewhat inferior results compared to the trivial+PDE preprocessing. Also, adding the full initial stress tensor and initial displacements to the PDE preprocessing did not improve over the presented PDE preprocessing.
We also experimented with different network architectures. In addition to the UNet, we considered a ResNet [54] of depth 10 with 5-layer CNN blocks. However, despite the UNet requiring less VRAM and computation time, it reliably outperformed the ResNet.

Conclusion.
We developed a physics-inspired deep learning approach to topology optimization. Afterward, we used a large-scale comparison to prove the effectiveness of our approach's two critical components -on the one hand, PDE-based preprocessing and, on the other, mirror and rotation equivariant architectures. These two key elements not only drastically improve the training's sample efficiency but also the predictions' physical correctness.
In practice, the existence of appropriate large-scale training data is costly or simply not realistic. Therefore, it is our conviction that sample efficiency is a crucial component of any future deep learning approach to topology optimization.
To promote future progress and comparability of the field, we also publish two 3D topology optimization datasets with a total count of almost 10 000 problem-ground truth pairs.       Our UNet architecture is inspired by [13] and based on the implementation in [35]. Like in [13], each layer of our encoder contains two 3×3×3 convolutions, each followed by batch normalization [24] and a rectified linear unit (ReLU) activation function, and then a max pooling operation. For the sphere dataset, we choose a pooling kernel size of 2 × 2 × 2; for the disc dataset, we choose 2 × 2 × 1 to account for the lower number of voxels in the z-direction. In the decoding path, each layer consists of a nearest neighbor upsampling step, followed by two 3 × 3 × 3 convolutions with batch normalization and ReLU activation functions. In the last layer, a 1×1×1 convolution with a sigmoid activation function reduces the number of output channels to 1 and acts as a binary voxel-wise classifier. We use skip connections to pass information from layers of equal resolution in the encoding path to the decoder and apply a dropout rate of 0.3 to reduce overfitting. We use five encoding and four decoding blocks. Each encoding block doubles the number of channels, starting with 10 output channels for the first block. We choose the padding, so the image resolution stays unaffected by the convolutions. Based on memory availability, we set the batch size to 128 when training on the disc dataset and 28 when training on the sphere dataset. All models have a total parameter count of approximately 1.6M. For the number of training epochs and the approximate training times for different models trained on the largest training subset, see Table 8. Appendix B. Convex hull preprocessing. For an alternative preprocessing strategy we construct a binary mask to extract all voxels with either Dirichlet boundary or load cases. The convex hull of this mask is then given as the set of voxels included in the smallest convex polygon that surrounds all these voxels. In case of one load point and one Dirichlet voxel this results in the simplest possible density distribution that leads to a connected mechanical structure. This convex hull density distribution is then used as a 1-channel input to the neural network.