Real-time design of architectural structures with differentiable mechanics and neural networks

1. In Figure 1, it seems that 3 figures are simultaneously provided to the encode, but the caption indicates only one shape is provided. Please take care to avoid confusion.

Thank you for noting this important detail. We have updated Figure 1 to match the caption, and to indicate clearly that one shape is fed to the encoder.

2. Describe the dimension of the matrix K.

Certainly! We have explicitly indicated the dimension of matrix $\mathbf{K}$ in Section 2.1 and in Appendix D, where we detail the math behind the mechanical simulator.

3. In Figure 3, it is hard to see how the x vectors are being updated.

We have retouched Figure 3 to improve legibility. We increased the size of the zoomed-in window (i.e., the callouts) in each step of the optimization, decreased the opacity of the background images in the callouts, and thickened the orange arrows that indicate the updates of $\mathbf{x}_i$.

6. Could you explain why the masonry shell is modeled as a pin-jointed structure but is then physically built as a shell (and how the two structure types are related)?

Excellent question! That might not be intuitive at first, but there is a substantial body of research in mechanics behind this modeling approach, e.g. [1, 2, 3]. This body of research revolves around two main concepts: limit state analysis and the Heyman safe theorem.

The behavior of masonry structures is governed by force equilibrium rather than material strength [1, 2]. This makes reasoning about the masonry behavior compatible with that of a pin-jointed bar system. The bar system hence represents a virtual network of compression forces flowing in the shell. This force network corresponds to one force state (among potentially many) in the masonry structure. If the force state satisfies the Heyman assumptions [2]; then the masonry structure is safe under the applied loads, and the solution corresponds to a lower bound of the collapse state (limit analysis).

The Heyman assumptions are [2, 3]:

• (i) Masonry’s compressive strength is infinite,
• (ii) masonry’s tensile strength is considered null, and
• (iii) no sliding occurs between the masonry blocks of the structure.

We have added a sentence to the paper at the start of Section 4.1 indicating these assumptions to the reader.

[1] S. Huerta, “The Analysis of Masonry Architecture: A Historical Approach: To the memory of Professor Henry J. Cowan,” Architectural Science Review, vol. 51, no. 4, pp. 297–328, Dec. 2008, doi: 10.3763/asre.2008.5136. [2] J. Heyman, “The stone skeleton,” International Journal of Solids and Structures, vol. 2, no. 2, pp. 249–279, Apr. 1966, doi: 10.1016/0020-7683(66)90018-7. [3] R. Maia Avelino, A. Iannuzzo, T. Van Mele, and P. Block, “Assessing the safety of vaulted masonry structures using thrust network analysis,” Computers & Structures, vol. 257, p. 106647, Dec. 2021, doi: 10.1016/j.compstruc.2021.106647.

7. In practice, how does one come up with a design $\hat{\mathbf{X}}$?

In practice, a design $\hat{\mathbf{X}}$ comes from an architect who proposes the global visual intent of a structure. Then, a structural engineer must counterpropose the architect’s proposal with $\mathbf{X}$, a mechanically efficient version of $\hat{\mathbf{X}}$. Broadly speaking, the level of happiness of the architect and the engineer is directly proportional to the distance between $\hat{\mathbf{X}}$ and $\mathbf{X}$. The closer the shapes are, the merrier. We envision the end users of our model to be both architects and engineers, facilitating the dialogue between both parties in a CAD environment at an accelerated pace (Appendix A).

8. How does the simulator only require p and b for shape prediction? It seems like the simulator is missing information to build a shell with N pins.

The information to calculate geometry for a system with $N$ pins is embedded in the structure of the stiffness matrix $\mathbf{K} \in \mathbb{R}^{N \times N}$, which is created by our mechanical simulator. The stiffness matrix is assembled by the simulator based on the connectivity between the $N$ pins and the $M$ bars. Since we keep this discretization and the connectivity constant in each task, the assembly information remains fixed across our experiments ($\mathbf{q}$, $\mathbf{p}$, and $\mathbf{b}$, in contrast, change values per every shape input to our model).

The initial submission included the equations to create the stiffness matrix, and we have now incorporated a more explicit description of the stiffness matrix assembly in Appendix D of our revised submission.

1. Are the authors going to be releasing the code?

Yes, we plan to release the code with the camera-ready version of the paper.

2. Does the chosen bar stiffness affect the loading on the design since a stiffer bar might lead to a heavier design?

Not necessarily. While stiffness and mass are related concepts, the loading is independent of the stiffnesses in our experiments. The reason is that stiffness is a metric that quantifies the ratio between the force passing through a given bar ($f$) and the length of the bar ($l$). See Equation 8 in our updated submission. Therefore, a high stiffness means that either the force in a bar is high (potentially requiring a thicker bar if the material the bar is made of is weak), or simply that the bar is very short (small $l$).

To avoid ambiguity, we define the magnitude of the loads upfront in our experiments. This is standard in engineering practice. In the shell task, for example, we set the magnitude to an area load of 0.5 force/area units (e.g., 0.5 kN / m2). This magnitude is commensurate with the loads in the materialized design generated by our model as a (tabletop) masonry vault of uniform thickness.

3. Could you provide further clarification on how the stress direction is determined a-priori?

We set the stress direction by setting the signs of the stiffness values $\mathbf{q}$ that we input to the mechanical simulator. We define the signs of $\mathbf{q}$ in the last layer of our encoder with the direction vector $\mathbf{s}$ (see Equation 3): a negative entry $s$ prescribes compressive stress, and a positive factor, tensile stress. We have added this clarification –-which was missing– to our updated submission, in L195-L196.

It is worth emphasizing that we can set the signs in the stress response of a structure a priori because our choice of specialized mechanical simulator allows us to do so for free (Appendix D). Defining the signs of the stress response a priori would be more significantly challenging with other traditional but still differentiable mechanical simulators, e.g. references [1] and [2], and would potentially involve equality constraints on the outputs of these other simulations. We mention this distinction in our description of the mechanical decoder in Section 2.3 in the updated submission.

[1] G. Wu, “A framework for structural shape optimization based on automatic differentiation, the adjoint method and accelerated linear algebra,” Struct Multidisc Optim, vol. 66, no. 7, p. 151, Jun. 2023, doi: 10.1007/s00158-023-03601-0.

[2] T. Xue et al., “JAX-FEM: A differentiable GPU-accelerated 3D finite element solver for automatic inverse design and mechanistic data science,” Computer Physics Communications, vol. 291, p. 108802, Oct. 2023, doi: 10.1016/j.cpc.2023.108802.

4. When training the neural network component, how is this loss accounting for permutation invariance with locations, or is there some implicit consistency in the ordering of the pin coordinates.

Our shape loss does not encode any permutation invariance, it only focuses on pointwise distance minimization between nodes (L159 - L161). To maintain consistency in our predictions, we keep the same order of the vertices in the rows in $\mathbf{X}$ in each task, as the reviewer correctly pointed out.

6. What mechanism is in place to prevent the model providing excessively large stiffness values for all the bars?

Our model has no explicit mechanism to limit the entries of $\mathbf{q}$. In the second task, we regularize training by adding a term to the loss function (Equation 6). Nevertheless, this regularization reduces the variance of $\mathbf{q}$ without imposing a specific upper bound. In direct optimization, we do enforce box constraints on stiffness values, capping $\mathbf{q}$ at 20. This information has been included in Appendix E of the revised manuscript.

1. Even for the shells, if (number of bars) and (the number of nodes) change, the model may have to be retrained. It's understandable that this does not reduce the practical value of the proposed method.

Our model would have to be retrained if the number of bars or nodes changes. This is because we currently use an MLP as the encoder, which fixes the dimensions of the input shape matrix X. However, we are not restricted to that encoder architecture since our model construction is general (L187-188).

We state this limitation in Section 6.1 (L533-534) and point to looking at other encoder architectures in future works to generalize across multiple bar connectivities without retraining (L534-535). However, as the reviewer suggests, the limitation does not reduce the value of our work.

To be more forward with the reader about our model’s limitations, we have reworded a pair of sentences in Section 6.1 in the updated version of our manuscript. We would also like to point out that training our model is reasonably fast (for example, 1.20 minutes in the shell task in Section 4.1), and that retraining would not be unreasonable if the connectivity changes to a moderate extent.

2. The runtime measurements have such large variances that they are not meaningful. The author(s) should consider changing the evaluation platform (macOS has many background processes and is not a good platform for benchmarking), and deploy better runtime measurement libraries.

Thank you for catching this problem. We build on our response to Reviewer NVMt to a similar question. The large variances in both tasks arose from a bug in our timing protocol. This bug was independent of the operating system.

Despite having pre-compiled the models, the just-in-time (JIT) compilation of JAX was re-triggered when we ran inference on the first shape in the test batch due to an inadvertent change in the batch dimensions. We were thus including the re-compilation time (in the 300-400 ms range) in our statistics, which drove the variances up.

We have fixed this bug and re-timed test inference for our model and the neural baselines. We have updated tables 1 and 2 with the new statistics in both tasks. Regarding tooling, now we have put checks in our code to ensure that JIT re-compilation time does not leak into our statistics again (e.g., by using tools like equinox.debug.assert_max_traces and setting the max number of traces to 1).

The variance is now in the submillisecond levels. After discarding the JIT compilation time, the NN and PINN inference times decreased by one order of magnitude, making for stronger baselines.

3. In line 200 in Eq. (3), the fact that $\tau$ can be used to specify minimal stiffness is very interesting. But there are no more follow-up discussion on this topic. Provide more discussion on the usage of $\tau$.

Thank you for the suggestion. While we briefly mentioned in our original submission (L203-202) that setting $\tau>0$ is equivalent to setting a lower box constraint on the stiffness values like in traditional optimization, we agree that the paper benefits from more discussion. To this end, we have extended Section 2.3 in our updated submission and mention that lower bounding the stiffness space with $\tau$ is beneficial for two reasons. First, it is beneficial to guide learning towards particular solutions since the map from $\mathbf{q}$ to $\mathbf{X}$ is non-unique due to mechanical indeterminacy [1, 2]. That is, multiple distributions of $\mathbf{q}$ can map to the same $\mathbf{X}$. Second, it is beneficial to provide numerical stability when learning to amortize the design of more complex systems that mix tensile and compressive bar forces, as we do in Task 2 in the paper. We also now state the value of tau we use in the shells task.

[1] T. Van Mele, L. Lachauer, M. Rippmann, and P. Block, “Geometry-Based Understanding of Structures,” Journal of the International Association for Shell and Spatial Structures, vol. 53, no. 4, pp. 285–295, Dec. 2012.

[2] S. Pellegrino and C. R. Calladine, “Matrix analysis of statically and kinematically indeterminate frameworks,” International Journal of Solids and Structures, vol. 22, no. 4, pp. 409–428, Jan. 1986, doi: 10.1016/0020-7683(86)90014-4.

1. Reporting training times would enable a better assessment of the method’s real-world viability; if training significantly exceeds the duration of several optimizations, direct optimization might be preferable.

We reported the training times of our model in Table 3 in the Appendix in our original submission (now Table 4). The training time of our model is in the order of a few minutes per task on CPU (140 seconds for the shells task and 810 for the towers task). Regarding viability, for example, the training time of our model is equivalent to solving 25 direct optimization problems in sequence in the shells task (taking 5.81 seconds as the average for optimization, Table 1). Once trained, our model allows for many more shell design iterations than 25.

2. Similarly, providing training time metrics for the other baselines would help clarify trade-offs and be useful in scenarios where slight accuracy losses are acceptable for performance gains.

We also provided training times for the neural baselines (NN and PINN) in Table 3 in the Appendix of our first submission (now Table 4). The training times of the baselines are in the order of minutes on a CPU but are shorter than our model because they do not run the mechanical simulation during training. We note that the difference in training times varies with the task complexity: in the shell task, the difference between the baselines and our model is only 1 minute, and in the cable net task (a more challenging task), that number rises to 11 minutes. While the NN and the PINN are faster to train, neither of them drives the physics loss to zero, which is one of the key drivers in our paper. Although our model is more expensive to train than the purely neural baselines, it satisfies the physics objective by construction and compensates for the training time by generating fast physics-constrained design iterations in inference in real time.

3. From Fig.9 , the optimization initialized with the proposed method converges quickly as compared to the other initializations. This leads to an interesting question of how the NN initialized optimizations would perform. This approach would have the benefit of outputting a guaranteed local minima while avoiding the additional physics overhead during training. Including these results would enhance the paper's contribution to the community.

Great idea! We performed an additional optimization experiment in the towers task, utilizing the PINN predictions as the initialization. We show the results in the main text and the updated version of Figure 9. Just like our model, the PINN initializations converge faster than expert initialization and arrive at better-performing designs. Nevertheless, the PINN initialization is slower to converge than our model and only reaches our model’s shape loss at the end of the convergence curve.