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

## Abstract

Abstract Designing mechanically efficient geometry for architectural structures like shells, towers, and bridges is an expensive iterative process.
Existing techniques for solving such inverse mechanical problems rely on traditional direct optimization methods, which are slow and computationally expensive, limiting iteration speed and design exploration.
Neural networks would seem to offer a solution, via data-driven amortized optimization for specific design tasks, but they often require extensive fine-tuning and cannot ensure that important design criteria, such as mechanical integrity, are met.
In this work, we combine neural networks with a differentiable mechanics simulator to develop a model that accelerates the solution of shape approximation problems for architectural structures modeled as bar systems.
As a result, our model offers explicit guarantees to satisfy mechanical constraints while generating designs that match target geometries.
We validate our model in two tasks, the design of masonry shells and cable-net towers.
Our model achieves better accuracy and generalization than fully neural alternatives, and comparable accuracy to direct optimization but in real time, enabling fast and sound design exploration.
We further demonstrate the real-world potential of our trained model by deploying it in 3D modeling software and by fabricating a physical prototype.
Our work opens up new opportunities for accelerated physical design enhanced by neural networks for the built environment.

## 1 Introduction

Mechanical efficiency is required for architectural structures to span hundreds of meters under extreme loads safely while minimizing mass and embodied carbon.
Additionally, shells, towers, and bridges—examples of such systems—must comply with geometric constraints arising from architecture and fabrication requirements to become feasible structures in the built environment.
Designing shapes for such long-span structures, which must fulfill mechanical efficiency and geometric constraints, is a complex task requiring substantial domain expertise and human effort.
Our goal is to use machine learning to accelerate this challenging task without compromising safety-critical aspects of the design.

One way to approach this problem is to start from the mechanical standpoint, employing a specialized mechanical model that directly computes efficient geometry for structures modeled as a bar systems (Bletzinger and Ramm, ; Bletzinger et al., ).
Unlike standard, finite-element-based mechanical analysis, where one first defines the structure’s geometry and then obtains its internal stresses, these specialized models – known as *form-finding methods* in structural engineering (Veenendaal and
Block, ; Adriaenssens et al., )– reverse the relationship between geometry and stress to produce mechanically efficient shapes in a forward solve (Shin et al., ).
As a result, these methods have been successfully applied to design landmark structures with thickness-to-span ratios up to 1:70 (less than that of an eggshell), across a wide palette of materials, including stainless steel (Schlaich, ), reinforced concrete (Isler, ), and stone (Block et al., ).

While utilizing a specialized model can lead to efficient designs, it is difficult to guide solutions toward particular geometries, as the designer only has explicit control of the mechanical behavior and not the shape.
To solve this inverse problem, form-finding is complemented with optimization algorithms to find stress states that satisfy geometric goals (Panozzo et al., ; Maia Avelino et al., ).
For a designer, however, it is cumbersome to perform a time-consuming and computationally intensive optimization when exploring shapes.
Practical design requires the evaluation of multiple target shapes to align with geometric desiderata, multiplying computational effort as each shape requires its own optimization.
Therefore, using specialized mechanical models and optimization together is effective but inefficient in practice, where shape variety and real-time feedback are essential.

Figure: Figure 1: Architecture of our model to amortize the generation of mechanically efficient geometry. Given an input target shape $\hat{\mathbf{X}}$, a neural network maps it to a stiffnesses space $\mathbf{q}$. The stiffnesses, in tandem with boundary conditions $\mathbf{b}$, are then decoded by a mechanical simulator into a physics-constrained shape $\mathbf{X}$ that approximates the input. The generated shape can then be used as the backbone geometry to design mechanically efficient structures like masonry shells.
Refer to caption: extracted/5899952/figures/neural_fofin_pipeline_5.png

Neural networks (NNs) have shown the potential to accelerate physical design with data-driven surrogate models that amortize inverse problems for more responsive tools.
Recent applications include fluid-structure control (Allen et al., ) and additive manufacturing (Sun et al., ); and specific to mechanical optimization, the design of truss lattices (Bastek et al., ), tall buildings (Chang and Cheng, ), reticulated shells (Tam et al., ) and cable-nets (Mai et al., ).
However, these purely data-driven approaches require the representation of both the inverse problem and the underlying physics.
Even with physics-informed neural networks (PINNs) (Raissi et al., ; Karniadakis et al., ) that have specialized architectures (Bastek and Kochmann, ; Lu et al., ), or that are trained with sophisticated loss balancing schemes (Bischof and Kraus, ; Wang et al., ), there is no guarantee of mechanical integrity in their predictions.
Assurance of mechanical integrity is a foundational tenet in structural design, where poor neural predictions might lead to catastrophic collapse and the loss of human lives.
In contrast, mechanical simulators in structural engineering have been developed for decades and offer a principled and interpretable way to model the physics of long-span structures.
These models capture the physics by construction.
A hybrid solution seems ideal, in which neural network amortization is integrated with differentiable physics models (Belbute-Peres et al., ; Thuerey et al., ; Um et al., ; Oktay et al., ; Yang et al., ) to construct a rich class of machine learning models that shift the current paradigm from a physics-informed to a physics-in-the-loop approach in safety-critical applications.

In this paper, we develop a neural surrogate model that couples a neural network with a differentiable mechanics simulator to enable the solution of shape approximation problems for architectural structures in real time (Fig. 1).
The coupled model offers advantages over direct gradient-based optimization and current fully neural alternatives for interactive mechanical design.
We evaluate our method in two design problems of increasing complexity: masonry shells and cable-net towers.
Our contributions are threefold.
First, we demonstrate that our model generates mechanically sound predictions at higher accuracy than NNs and PINNs of similar architecture.
The model exhibits better generalization performance than an equivalent PINN in the masonry shells task.
Second, we show that our model reaches comparable accuracy to optimization, but our model is on average three orders of magnitude faster in both design problems.
In the cable-net task, our model provides robust initialization for direct optimization, outperforming the designs generated by optimization initialized with human domain expertise.
Third, we showcase the application of a maturing machine learning technique (i.e., coupling learnable and analytical components in the same architecture) to a new, high-impact domain for physical design (i.e., architectural structures at the meter scale).
To illustrate its practical impact, we deploy our trained model in a 3D modeling program to solve a shell design problem and then fabricate a physical prototype of the predicted geometry.

## 2 Physics-Constrained Neural Form Discovery

Figure: Figure 2: A bar system. In the callout, the stiffness component $[\mathbf{K}(\mathbf{q})]_{ij}$ is indicated at a bar connecting nodes with positions $\mathbf{x}_{i}$ and $\mathbf{x}_{j}$. A load $\mathbf{p}_{i}$ is applied at $\mathbf{x}_{i}$. In this system, the anchor nodes on the perimeter are fixed.
Refer to caption: extracted/5899952/figures/fdm_setup.png

Our goal is to generate mechanically efficient shapes for architectural structures that approximate target geometries in real time while maintaining mechanical integrity.
The challenge is that, because mechanical integrity is paramount in architectural structures, it is necessary to reason about designs from the point of view of force balance; but the resulting geometries are a nontrivial function of their mechanical behavior.
Thus we seek to use machine learning to efficiently invert this function to generate target designs without compromising their integrity.

### 2.1 Form-finding efficient geometry

We focus on structures modeled as pin-jointed bar systems of $N$ nodes connected by $M$ bars (Figure 2).
Each node experiences an external load vector (e.g., self-weight and wind load) and some nodes are constrained to fixed positions (e.g., terrain and anchors).
These are the structure’s boundary conditions $\mathbf{b}\in\mathbb{R}^{L}$.
After picking bar stiffnesses $\mathbf{q}\in\mathbb{R}^{M}$, the goal of form-finding is to identify positions of the nodes $\mathbf{X}=(\mathbf{x}_{1},\ldots,\mathbf{x}_{N})\in\mathbb{R}^{N\times 3}$ such that there is no net residual force on the structure.

The force density method (Schek, ) is a canonical form-finding method that works on structures modeled as bar systems.
This method generates mechanically efficient shapes by mapping the bar stiffnesses $\mathbf{q}$ to an overall shape represented by the node positions $\mathbf{X}$ in a state of equilibrium.
Appendix D offers a more thorough description.
An efficient structure is one whose loaded configuration is bending and torsion-free.
This configuration minimizes the structure’s total strain energy by reducing the contribution of such components, letting a structure sustain applied loads mainly under tensile and compressive axial stresses.
To arrive at equilibrium and satisfy physical laws, the residual force vector $\mathbf{r}_{i}\in\mathbb{R}^{3}$ for a free node with position $\mathbf{x}_{i}$ must be zero.
The function $\mathbf{r}_{i}(\mathbf{X};\mathbf{q})$ quantifies the difference between the load $\mathbf{p}_{i}\in\mathbb{R}^{3}$ applied to node $i$ and the sum of the internal stresses of the bars incident to the node, for given node positions, bar stiffnesses and boundary conditions:

$$ $\mathbf{r}_{i}(\mathbf{X};\mathbf{q}):=\sum_{j\in\mathcal{N}(i)}[\mathbf{K}( \mathbf{q})]_{ij}(\mathbf{x}_{i}-\mathbf{x}_{j})-\mathbf{p}_{i}(\mathbf{x}_{i})$ (1) $$

where $\mathcal{N}(i)$ are the neighbors of node $i$, and $\mathbf{K}(\mathbf{q})$ is the stiffness matrix as a function of $\mathbf{q}$.
Positive and negative stiffnesses correspond to tension and compression stresses, respectively.
The restriction that there is no net residual force can be framed as a linear system in which all of the $\mathbf{r}_{i}(\mathbf{X};\mathbf{q})=\mathbf{0}$.
Although it can be solved using conventional direct methods, the resulting map from the stiffnesses $\mathbf{q}$ to the positions, which we denote $\mathbf{X}(\mathbf{q})$ is implicit and nonlinear, so, difficult to reason about directly.

### 2.2 Direct optimization for target shapes

Although form-finding methods have the appealing property that they guarantee mechanical efficiency, they do not allow a designer to directly target particular geometries.
Moreover, not all geometries are even compatible with mechanical efficiency.
If a designer has a target shape $\hat{\mathbf{X}}$, they wish to solve the following optimization problem with respect to $\mathbf{q}$ to approximate $\hat{\mathbf{X}}$ with $\mathbf{X}$:

$$ $\displaystyle\mathbf{q}^{\star}$ $\displaystyle=\operatorname*{arg\,min}_{\mathbf{q}\in\mathbb{R}^{M}}\mathcal{L }_{\text{shape}}(\mathbf{q})$ where $\displaystyle\mathcal{L}_{\text{shape}}(\mathbf{q})$ $\displaystyle:=\sum_{i=1}^{N}\sum_{d=1}^{3}|[\mathbf{X}(\mathbf{q})]_{i,d}-[ \hat{\mathbf{X}}]_{i,d}|^{p}\,,$ (2) $$

and $p>0$.
We call $\mathcal{L}_{\text{shape}}$ the *shape loss*.
This objective, which we refer to as *direct optimization*, tries to identify stiffnesses $\mathbf{q}$ which are close to the designer’s intent in an $\ell_{p}$ sense, while maintaining net zero force balance.
Note the optimization setup does not contain any information about the set of physically valid forms and it is simply driven by the minimization of the pointwise difference between shapes.
Conventionally, this nonlinear optimization problem needs to be solved numerically in the inner loop of a design process, but that is slow and computationally costly.

Figure: Figure 3: Inverse form-finding. To generate a mechanically congruent shape $\mathbf{X}(\mathbf{q})$ that approximates an arbitrary target $\hat{\mathbf{X}}$, traditional methods like direct optimization find bar stiffnesses $\mathbf{q}$ that minimize the shape loss $\mathcal{L}_{\text{shape}}$, but only after several iterations $t$. Our model amortizes this computationally taxing process during inference while guaranteeing physics, enabling real-time and sound design.
Refer to caption: x1.png

### 2.3 Amortized inverse form-finding

Rather than performing many costly optimizations within a design loop, we can use a machine learning model to amortize the solution over a family of target shapes $\hat{\mathcal{X}}$.
Specifically, we construct a neural network that takes as input the designer’s intent $\hat{\mathbf{X}}$, and outputs a set of stiffnesses $\mathbf{q}$ such that $\mathbf{X}(\mathbf{q})\approx\hat{\mathbf{X}}$.
As shown in Figure 1, our model architecture resembles an autoencoder.
First, a neural encoder $\mathcal{E}_{\phi}$ maps the target shape $\hat{\mathbf{X}}$ into a $\mathbf{q}$.
Then, a differentiable mechanical simulator implementing the force density method *decodes* the associated shape $\mathbf{X}(\mathbf{q})$, conditioned on the boundary restraints $\mathbf{b}$.
The key property of this construction is that the resulting shape is mechanically sound even if the neural network is inaccurate; the failure mode is not a lack of structural integrity, but a shape that does not match the targets very well.

##### Neural encoder

The encoder is a neural network with learnable parameters $\phi$ that ingests the targets and projects them into stiffness space, $\mathcal{E}_{\phi}:\mathbb{R}^{3N}\to\mathbb{R}^{M}$.
For simplicity, we cast our problem as a point-wise matching task and use a multilayer perceptron (MLP) as the encoder, although we are not restricted to that.
Regardless of the neural network specification, the output representation must be strictly positive.
This is necessary for compatibility with our mechanical simulator to avoid null stiffness values that bear limited physical meaning in our representation of an architectural structure and lead to singular stiffness matrices.
Here, we satisfy this requirement by applying a strictly positive nonlinearity to the last layer of the encoder.

One of the advantages of our mechanical simulator is that it enables us to prescribe tensile or compressive bar stresses *a priori*.
As a result, rather than making these stress directions a learnable feature, we build this bias into the encoder architecture by scaling the strictly positive embedding of the last layer by a stress direction vector $\mathbf{s}\in\mathbb{R}^{M}$.
The scaling factors $s_{i}\in\{-1,1\}$ indicate the direction of the internal axial stress of every bar.
Our encoder thus calculates a stiffness vector as:

$$ $\mathbf{q}=\mathbf{s}\odot(\sigma(\mathbf{h})+\tau)$ (3) $$

where $\sigma$ is the strictly positive nonlinearity, $\mathbf{h}$ is the encoder’s last layer embedding, ${\tau>0}$ is a fixed scalar shift that specifies a minimum absolute stiffness value for the entries of $\mathbf{q}$ (akin to box constraints in numerical optimization), and $\odot$ indicates element-wise product.

##### Training

To amortize the shape-matching problem, we look for model parameters $\phi^{\star}$ that minimize the expected value of the shape loss over a family of target shapes $\hat{\mathcal{X}}$:

$$ $\phi^{\star}=\operatorname*{arg\,min}_{\phi}\mathbb{E}_{\hat{\mathbf{X}}\sim \hat{\mathcal{X}}}\left[\sum_{i=1}^{N}\sum_{d=1}^{3}|[\mathbf{X}(\mathcal{E}_{ \phi}(\hat{\mathbf{X}}))]_{i,d}-[\hat{\mathbf{X}}]_{i,d}|^{p}\right]\,.$ (4) $$

We train our model via first-order stochastic gradient descent, averaging the loss values over batches at each training step.
Appendix E provides training specifications.
We generate training data by sampling batches of target shapes from a task-specific design space parametrized by a probability distribution (see Section 4).
Our model can be trained end-to-end, coupling the neural network and the mechanical simulator, because the encoder and decoder are both implemented in a differentiable programming environment (Bradbury et al., ).
As a result, reverse-mode automatic differentiation can seamlessly backpropagate the physics-based gradients that tune the neural network parameters.

## 3 Evaluation

We evaluate model performance by measuring the inference wall time of the trained model in addition to the value of the shape loss $\mathcal{L}_{\text{shape}}$ over a test batch of shapes of size $B$.
We compare the performance of our model to three other baselines: a fully neural approach (NN), a fully neural model augmented with a physics-informed loss (PINN), and traditional direct optimization.
The fully neural approach replaces the differentiable simulator in our model with a learnable decoder mirroring the encoder’s architecture, with the inclusion of the boundary conditions $\mathbf{b}$ as inputs, as is the case with the differentiable simulator.
The fully neural model is then trained to minimize the shape loss $\mathcal{L}_{\text{shape}}$ without any additional regularization, highlighting that an information bottleneck is insufficient to automatically guarantee physically plausible designs.
The second baseline extends the fully neural approach by adding an explicit physics loss:

$$ $\mathcal{L}_{\text{physics}}=\mathbb{E}_{\hat{\mathbf{X}}\sim\hat{\mathcal{X}} }\left[\sum_{i=1}^{N}||\mathbf{r}_{i}(\mathbf{X}(\mathcal{E}_{\phi}(\hat{ \mathbf{X}})))||_{2}\right]$ (5) $$

The physics loss $\mathcal{L}_{\text{physics}}$ is the governing equation of our problem, and it measures the $N$ residual forces $\mathbf{r}_{i}$ in the nodes.
For a shape to be mechanically sound in our setup, $\mathcal{L}_{\text{physics}}$ must be zero within numerical precision (i.e., $1\times 10^{-12}$) as per our physics constraint in Equation 1.
We reason that the additional term should provide a training signal to the encoder and decoder such that they learn how to solve the shape-matching tasks and the physics concurrently.
The third baseline takes advantage of the differentiable physics simulator and directly optimizes the parameter space $\mathbf{q}$ input to the decoder to minimize the shape loss via deterministic gradient-based optimization on a per-shape basis.

## 4 Experiments

We test our model to amortize shape approximation tasks for masonry shells and cable-net towers.
These two tasks represent a broad class of structural typologies dealt with by designers in practice.

Figure: (a) Train curves for $\mathcal{L}_{\text{shape}}$ and $\mathcal{L}_{\text{physics}}$
Refer to caption: x2.png

### 4.1 Masonry shells

Our first experiment identifies suitable bar stiffness values for unreinforced masonry shells.
Masonry shells sustain external loads with span-to-thickness ratios as low as 1:40 despite being built from materials that are strong in compression and weak in other loading conditions (Block et al., ).
Structural forms maximizing internal compressive axial stresses enable this efficient behavior.

Figure: (a) NN
Refer to caption: x5.png

**Table 1: Model evaluation on masonry shells task. We report mean loss values and the standard deviation per shape on the test set of 100 target shapes and inference run time.**
|  | NN | PINN | Optimization | Ours |
| --- | --- | --- | --- | --- |
| $\mathcal{L}_{\text{shape}}\,\downarrow$ | $1.4~{}\pm~{}0.3$ | $3.1~{}\pm~{}1.4$ | $0.8~{}\pm~{}1.3$ | $3.0~{}\pm~{}2.0$ |
| $\mathcal{L}_{\text{physics}}\,\downarrow$ | $97.3~{}\pm~{}35.6$ | $0.7~{}\pm~{}0.3$ | $0.0~{}\pm~{}0.0$ | $0.0~{}\pm~{}0.0$ |
| Time [ms] $\downarrow$ | $3.3~{}\pm~{}30.1$ | $3.3~{}\pm~{}30.1$ | $5810.1~{}\pm~{}1870.3$ | $4.1~{}\pm~{}34.7$ |

An expressive class of masonry shells can be constructed by surfaces parameterized by a Bezier patch.
The shape of the patch is in turn described by the position of a grid of $C$ control points $\mathbf{c}$ in Cartesian space (see Appendix F).
For this task, we restrict the space of target shapes to a square grid of width $w=10$ and $C=16$ control points.
In particular, we consider doubly-symmetric shapes of constant discretization, where $N=100$ and $M=180$.
We apply a constant area load of $0.5$ per unit area, representing the self-weight of the shell.
The nodes on the perimeter are anchored.

To solve this task, we look for bar stiffnesses that yield shapes that best fit the target geometries, and whose internal axial stresses are compressive ($\mathbf{q}<0$).
We satisfy the compression-only requirement by construction with our encoder by setting $\mathbf{s}=-\mathbf{1}$.
We use $p=1$.
Figure 5(a) illustrates the stochastic loss curves during training for our model and the two neural baselines (NN and PINN).
The fully neural approaches achieve a low shape loss but are unable to converge w.r.t the physics loss.
These approaches fail to learn a meaningful intermediary representation $\mathbf{q}$ that satisfies the task physics, unlike our physics-in-the-loop model where this requirement is satisfied by construction.
In Figure 8, we plot two representative target shapes for each model, the predicted stiffnesses, the resulting geometries, and the residual forces.

The trends we observe in training are repeated during inference.
Table 1 reports the average test loss and the standard deviation for one design in a test group of 100 different structures.
Direct optimization achieves the lowest shape loss.
The baselines and our model all offer significant speedup w.r.t optimization.
While the NN and the PINN generate accurate shape approximations, the designs predicted by these fully neural baselines are mechanically unfeasible because the residual forces fail to vanish.
The magnitude of the physics loss indicates that the structure is missing balancing forces to achieve equilibrium.
If we prescribe typical values for a masonry shell constructed out of bricks with an average area of $1$m^2 and thickness of $0.05$ m (a $1:20$ aspect ratio), density of approximately $2000$ kg/m^3 then a total residual force of 1 kN, representative of a unit value of $\mathcal{L}_{\text{physics}}$ in Table 1, would destabilize the structure with an acceleration of $10$ m/s^2 in multiple directions.
Consequently, the NN and the PINN shapes are unstable, and building masonry shells guided by these predictions can lead to collapse.
In contrast, our model and direct optimization satisfy the physical constraints a priori, but our model generates accurate predictions up to three orders of magnitude faster, and offers significant speedup compared to optimization for geometries with an equivalent shape loss value (Appendix B).

Figure: $\mathcal{L}_{\text{shape}}=16.6$
Refer to caption: x14.png

Next, we investigate the out-of-distribution generalization of the PINN and our model, a desirable property to build robust neural surrogates for physical design.
To this end, we first generate a set of 100 asymmetric Bezier surfaces.
We create this new set by sampling control points from the same design space and discretization as in the doubly-symmetric case (Appendix F).
Then, we produce input data by interpolating between the shapes in the asymmetric and symmetric sets.
We evaluate the shape loss of the model predictions at increasing interpolation factors $\delta$, where $\delta=0$ denotes double symmetry and $\delta=1$ indicates full asymmetry.
We do an equivalent study for $\mathcal{L}_{\text{physics}}$ in Appendix B.
Our model consistently possesses better out-of-distribution performance.
Figure 5(b) demonstrates that the loss of our model predictions decays at least at half of the PINN’s rate, and with a lower spread, as we increase the data asymmetry.
At $\delta=1$, the loss of our model is $2.5$ times lower.
The performance disparity at this point between our model and the PINN is evidenced by the example in the third column of Figure 19b.
The PINN prediction is not only a worse fit to the target, but also the residuals are fifteen times higher than at $\delta=0$, demonstrating that physics-in-the-loop models offer enhanced generalization over physics-informed approaches.

To further demonstrate the applicability of our model under real-world conditions, Figure 15 in Appendix A shows screen captures of our trained model in action, assisting a designer in exploring different shapes for a masonry shell in an industry-grade 3D modeling software.
We also validate our model’s generalization and transfer to physical applications by fabricating a tabletop masonry shell based on one of our model predictions on asymmetric targets displayed in Figure 10(a).
After tessellating the predicted shape, we manufacture the individual bricks and assemble them.
The prototype is stable and can resist gravity without glue or mechanical connectors, indicating that the structure can work predominantly under internal compressive forces as required by masonry shells, and proving the feasibility of our predictions under real-world constraints.
See Figure 10(b) for details.

### 4.2 Cable-net towers

Guided by the success of our initial experiment, we turn to a more complex problem: the design of cable-net towers.
Such systems are extremely lightweight, with low structural mass-to-occupied volume ratios, often being found as observation or cooling towers (Adriaenssens et al., ).
To explore the design space of these towers, we focus on optimizing the shape and orientation of horizontal compression rings encircled by vertically spanning tensile cable-nets.

Figure: (a) Task
Refer to caption: x15.png

Figure 12(a) gives an overview of the task setup.
The cable-net towers consist of two tensile nets that interface at a compressive ring at mid-height.
Each cable-net tower comprises $D=21$ rings with 16 points each, spaced at equal $h/D$ intervals over a total of height $h=10$.
The discretization of the structure is $N=335$ and $M=656$.
We parametrize the geometry of the bottom, middle, and top rings with an ellipse of radii $\alpha_{1}r$ and $\alpha_{2}r$ and rotation angle on the plane $\beta$.
We detail the data generation process in Appendix C.
The nodes on the top and bottom rings are anchors.
The shape approximation task here consists of finding valid cable-net shapes that match the target shape of the middle compression ring.
We are interested in geometries where the tension rings are planar.
The shape loss incorporates these requirements, measuring the squared $\ell_{2}$ norm (i.e. $p=2$) between predictions and targets.
We set the entries of $\mathbf{s}$ to $-1$ and $1$ to enforce the stress direction in the middle ring and the net bars, respectively.

Figure: Figure 13: Train loss curve of the cable-net task, with and without clipping.
Refer to caption: x20.png

The mechanical behavior of cable-net towers introduces modeling difficulties because the resulting problem to be solved comprises members that are either under tensile or compressive stresses.
This results in an ill-conditioned system due to the interaction between bars of different stress signs and can lead to singular systems at stress levels near zero (Cai et al., ).
While the ill-conditioning can be overcome in the forward pass, in the backward pass, it can hinder our model from learning meaningful representations as a result of poorly scaled gradients (Figure 13).
To alleviate these numerical instabilities, we clip the global gradient norm to $0.01$ and shift the outputs of our last layer to establish a lower bound of $\tau=1$ on the stiffness space.
Additionally, we add a regularization term to the total loss to steer our model towards a uniform stiffness distribution:

$$ $\mathcal{L}_{\text{reg}}=\text{Var}\,(\mathbf{q}_{\text{pos}})+\text{Var}\,( \mathbf{q}_{\text{neg}})$ (6) $$

The regularizer measures the variance, $\text{Var}(\mathbf{q})$, of the structure’s bar stiffnesses, for tensile $\mathbf{q}_{\text{pos}}$ and compressive $\mathbf{q}_{\text{neg}}$ values, over all the $B$ samples in a batch.
We scale $\mathcal{L}_{\text{reg}}$ by a constant factor $\lambda$.
We employ $\lambda=10$ to train our model and the baselines.

**Table 2: Model evaluation on the cable-net towers task. We report mean loss values and the standard deviation per shape on the test set of 100 towers and inference run time.**
|  | NN | PINN | Optimization | Ours |  |  |
| --- | --- | --- | --- | --- | --- | --- |
|  |  |  | Randomized | Expert | with Ours |  |
| $\mathcal{L}_{\text{shape}}\,\downarrow$ | $7.5~{}\pm~{}3.1$ | $7.4~{}\pm~{}3.4$ | $14.3~{}\pm~{}28.2$ | $0.3~{}\pm~{}0.3$ | $0.2~{}\pm~{}0.3$ | $0.6~{}\pm~{}0.7$ |
| $\mathcal{L}_{\text{physics}}\,\downarrow$ | $45.8~{}\pm~{}0.4$ | $1.2~{}\pm~{}0.1$ | $0.0~{}\pm~{}0.0$ | $0.0~{}\pm~{}0.0$ | $0.0~{}\pm~{}0.0$ | $0.0~{}\pm~{}0.0$ |
| Time [ms] $\downarrow$ | $1.6~{}\pm~{}11.4$ | $1.6~{}\pm~{}11.4$ | $1902.0~{}\pm~{}1297.3$ | $1876.2~{}\pm~{}821.5$ | $1430.1~{}\pm~{}838.9$ | $7.4~{}\pm~{}28.1$ |

Figures 12(b)-12(d) show an illustrative example of the predicted cable-net shapes.
In Appendix C, we show that our model predictions cover the task space satisfactorily, generating accurate and mechanically congruent cable-net shapes whose rings radii and in-plane rotation vary within a range of $r$ and $\pi/6$, respectively.
Table 2 demonstrates that our model generates shapes that match the targets with a 3 times tighter fit than the PINN and the NN models, although these models are faster than ours in generating shapes end-to-end.
Moreover, the neural baselines are unable to learn the cable-net physics because $\mathcal{L}_{\text{physics}}$ is nonzero within numerical precision, as we identified in the shell design task.

Next, we compare the two approaches that guarantee physics: ours and direct optimization.
Converging to a good local optimum with optimization is contingent on the choice of adequate initial stiffnesses $\mathbf{q}$.
In Figure 14 we thus analyze the expected convergence rate of direct optimization on the test set with three different initializations.
Random initial guesses lead to convergence to poor local optima as the shape loss is the highest among our experiments, highlighting the relevance of picking adequate $\mathbf{q}$ values to optimize the geometry of structures with complex mechanical behavior.
If the initial $\mathbf{q}$ values are handpicked by a human expert, optimization is more accurate than our model and the other baselines.
Optimization with expert initialization achieves an equivalent shape loss to our model’s at about $1/5$ of the total convergence time ($350$ ms) given by Table 2.
At this point, optimization has achieved most of the shape loss reductions w.r.t. the initial values, and further optimization only leads to diminishing returns on accuracy as illustrated by the decreasing slope of the convergence curve.
However, our trained model matches optimization in only 1 inference step that is on average three orders of magnitude faster and is free of potentially expensive human intervention—key attributes for fast and automated structural design tools.

Figure: Figure 14: Convergence rate of direct optimization with three distinct initializations.
Refer to caption: x21.png

Lastly, we investigate the effect of using our model predictions as the initial stiffness input to direct optimization to further refine the tower designs.
This combination results in the most accurate matches across all our cable-net tower experiments, converging faster and consistently achieving a lower shape loss than direct optimization with expertly initialized parameters (Figure 14).
That is, the neural model provides a better initial guess than the human expert for optimization in this task.
This finding suggests that neural networks and standard techniques are not mutually exclusive and that utilizing them in tandem is a promising direction to arrive at best-performing designs for architectural structures, and should be further explored in subsequent research.

## 5 Related works

##### Differentiable mechanical simulators

Machine learning and automatic differentiation have successfully obtained derivatives of complex forward mechanical models by either learning a differentiable surrogate with a neural network, or by implementing the analytical physics model in a differentiable programming environment.
For neural surrogates, Xue et al. () used an autoencoder to simulate 2D metamaterials, while Zheng et al. () learned the forward solution of 3D graphic statics on compression-dominant shells.
Differentiable simulator implementations for mechanics problems include the finite element method for 3D solids (Xue et al., ) and isogeometric analysis of 3D membranes (Oberbichler et al., ).
Specifically for form-finding, differentiable simulators exist for matrix structural analysis (Wu, ), and for combinatorial equilibrium modeling (Pastrana et al., ).
Our work utilizes one such form-finding simulator and connects it with a neural network to amortize an inverse problem.

##### Amortization models for mechanical design

Neural surrogates that amortize inverse problems have garnered attention for their ability to approximate nonlinear relationships in mechanical systems between target properties and input parameters.
At the centimeter scale, Bastek et al. () address the inversion of the structure-property map in truss metamaterials, enabling the discovery of optimal configurations.
In Oktay et al. (), amortized models were used to generate actuation policies to deform 2D cellular metamaterials to several target configurations.
Focusing on architectural structures at the meter scale, Hoyer et al. () proposed a neural basis for computing material distributions that minimize compliance for buildings in 2D; while Chang and Cheng () amortized with graph neural networks the design of the cross sections of the beams and columns of buildings in 3D against multiple load cases.
Favilli et al. () applied geometric deep learning and a differentiable simulator to gridshell structures for low-strain energy shapes.
Unlike their work, which trains a neural network for a single problem, we amortize over multiple inverse problems.
Tam et al. () also amortize a shape-matching problem with form-finding like in our work, but they do so using a fully neural approach that must capture both physics and solve the matching task simultaneously.

## 6 Conclusion

We presented a physics-in-the-loop neural model that expedites the solution of shape-matching problems to design mechanically efficient architectural structures.
By embedding prior physics knowledge in a neural network and training end-to-end, our model learns representations that solve a family of inverse problems with precision while enforcing mechanical constraints by construction; which is where current neural approaches fall short.
While we do not discard that physics-informed networks could push the physics loss closer to zero in the limit of extensive hyperparameter tuning and network size, the physics guarantees and the stronger generalization performance of our model set it apart as a more robust approach to support design in practical settings.
Our method generates mechanically sound geometries in milliseconds, with accuracy comparable to gradient-based optimization algorithms.
The speed and reliability of our model enable real-time design exploration of long-span structures, such as masonry shells and cable-net towers.

### 6.1 Limitations and future work

Although not a cure-all for mechanical design problems, our work evidences that domain expertise in structural engineering and machine learning is necessary to address numerical pathologies that can stem from infusing physics into neural networks (Wang et al., , ; Metz et al., ).
In our case, the simulator can yield a stiff ill-conditioned system which affects training in the presence of structures with complex stress distributions.
Tackling these instabilities requires knowledge of the physics of the problem being modeled, while also utilizing gradient stabilization techniques common in machine learning, such as gradient clipping.

Like past methods at this intersection, our model requires devising specialized parameter spaces for learning (Allen et al., ).
Consequently, one of our trained models may not necessarily generalize to different classes of architectural structures.
Another limitation is that our model is currently restricted to a single topology.
In the future, however, we expect that methods such as graph networks (Battaglia et al., ; Pfaff et al., ) will enable moving beyond a fixed bar connectivity, and developing more general models that accelerate the design of architectural structures across multiple classes.
We additionally note that the choice of bar stiffnesses for a given system is not unique and it is potentially appealing to present to the designer a diversity of possible solutions by reformulating our model in a variational setting (Kingma and Welling, ; Salamanca et al., ).
That is another exciting avenue for future research.

#### Acknowledgments

The authors thank Deniz Oktay for early discussions on this project.
This work has been supported by the U.S. National Science Foundation under grant OAC-2118201.

## Appendix A Real-time design in CAD software

We deploy our trained model in a 3D modeling software called Rhino3D [Robert McNeel & Associates, ] as illustrated in Figure 15.
Rhino3D supports traditional computer-aided design workflows and Grasshopper (Figure 15 top left), its visual programming extension, allows the creation of new software features via custom Python scripts.
We load our model as a Grasshopper plugin via Python and test it to support the real-time exploration of shapes for masonry shells.

We describe a design exploration session next.
A designer models a Bezier surface in Rhino3D by hand and imports it into Grasshopper.
Then, as the designer moves the control points of the Bezier, the geometry automatically updates.
Our model generates new compression-only shapes (see the bar systems rendered in Figure 15).
Note that as we detail in Section 4, we trained our model on fully symmetric surfaces but the adequate generalization it exhibits to asymmetric geometries makes it possible to support the designer during their exploratory session.

Figure: Figure 15: A demonstration of our trained model providing predictions of compression-only shapes in real time. The first screenshot shows the environment running our custom code. The screenshot in the second column and first row shows the initial input Bezier surface with control points. Every subsequent pair of screenshots shows the designer moving Bezier control points and our model reacting to the designer’s intent by approximating the shape with one that is compression-only.
Refer to caption: extracted/5899952/figures/rhino_bezier_4.png

## Appendix B Complementary studies for masonry shell design

Figure: (a)
Refer to caption: x22.png

### B.1 Comparison with direct optimization

Both our model and direct optimization satisfy the physical constraints a priori, but they generate predictions at drastically different speeds.
Since optimization is an iterative approach, we compare the evolution of the shape loss over run time to that of our method, which predicts the bar stiffnesses of a structure in one step.
Figure 16(a) shows that optimization reaches the best predictive performance among the baselines (Table 1), but only after convergence, which incurs run time over 5000 ms.
Even though optimization matches the shape loss of our model in $40\%$ of the expected convergence time, our approach offers a speedup of above $450\times$ as it generates designs of equivalent accuracy in only a few milliseconds.

### B.2 Generalization on the physics loss

In Section 4.1, we discussed the generalization capacity of our model and a PINN —the two approaches that are aware of the inverse problem physics— by observing the rate of change of the shape loss $\mathcal{L}_{\text{shape}}$ as we morph the geometry of the target shapes from doubly-symmetric to asymmetric.
Both models were trained solely on the former type.

We now quantify the impact of perturbing the test data distribution on the models’ capacity to preserve the problem’s physics.
Figure 16(b) depicts the evolution of the physics loss $\mathcal{L}_{\text{physics}}$ in logarithmic scale.
The performance of the PINN model rapidly deteriorates as the distribution changes from doubly-symmetric to asymmetric shapes, until the loss becomes one order of magnitude higher than at the start.
This finding illustrates that reasonable PINN performance in-distribution is not conducive to good generalization w.r.t. physics information out-of-distribution.
In contrast, our model satisfies the physics by construction, so $\mathcal{L}_{\text{physics}}$ remains constant and effectively zero within numerical precision regardless of the geometry of the target shapes.
Figure 19 depicts this trend with examples.

Figure: $\delta=0$
Refer to caption: extracted/5899952/figures/bezier_lerp_00_target_50_saddle_top.png

## Appendix C Shape exploration for cable-net tower design

### C.1 Data generation

To generate data for the cable-net task described in Section 4.2, we parametrize the geometry of the bottom, middle, and top rings with an ellipse of radii $\alpha_{1}r$ and $\alpha_{2}r$ and rotation angle $\beta$.
After setting the reference radius to $r=h/5$, we generate shape variations by sampling the scale factors $\alpha_{i}\in[1/2,3/2]$, and the angles $\beta\in[-\pi/12,\pi/12]$ from a uniform distribution.
See Figure 12(a) for a graphical description.

### C.2 Shape exploration

We carry out two experiments that modify the size and the rotation of the middle compressive ring in a cable-net tower to illustrate the geometric range of our model predictions.
In both experiments, we keep the scale of the radii of the top and bottom disks circular and equal by setting $\alpha_{1}=\alpha_{2}=3/4$.
In the first experiment, shown in Figure 21, we show our model predictions as we linearly scale the radius $\alpha$ of the middle circular ring from $\alpha=1/2$ to $\alpha=3/2$.
In Figure 23, the middle ring is elliptical ($\alpha_{1}=1/2$ and $\alpha_{2}=3/2$), and we twist it around the $z$ Cartesian axis in five steps between $-\pi/12$ and $\pi/12$.
The predictions in both tests match the target geometry.
Note that the distribution of the bar stresses $\mathbf{f}$ changes per structure due to the geometric changes on the middle ellipse.

Figure: $\alpha=1/2$
Refer to caption: x24.png

Figure: $\beta=-\pi/12$
Refer to caption: x29.png

## Appendix D Mechanical simulator

We employ the force density method (FDM) [Schek, ] as our mechanical simulator.
The FDM is a form-finding method that computes torsion and bending-free shapes on pin-jointed bar systems with $N$ nodes and $M$ bars, by mapping the bars’ stiffness vector $\mathbf{q}\in\mathbb{R}^{M}$ to node positions $\mathbf{X}\in\mathbb{R}^{N\times 3}$; subject to boundary conditions $\mathbf{b}\in\mathbb{R}^{L}$, $\mathcal{F}:\mathbb{R}^{M}\to\mathbb{R}^{3N}$.
In the FDM, a portion of nodes of size $N_{s}$ is fixed (i.e., anchors), and another of size $N_{u}$ is free to displace.
The size of the partitions is arbitrary as long as $N=N_{s}+N_{u}$.
For the geometry of a bar system to be in equilibrium and satisfy physical laws, the residual force vector $\mathbf{r}_{i}$ at each free node $i$ must be zero.
The residual quantifies the difference between a load applied to a free node $i$ and the sum of the $j$ incident bar axial stress vectors.
To satisfy this constraint, the FDM calculates the free node positions $\mathbf{X}_{u}\in\mathbb{R}^{N_{u}\times 3}$ after solving a linear system:

$$ $\mathbf{X}_{u}=\mathbf{K}(\mathbf{q})_{u}^{-1}\,[\mathbf{P}_{u}(\mathbf{X})- \mathbf{K}(\mathbf{q})_{s}\,\mathbf{X}_{s}(\mathbf{X})]$ (7) $$

where the stiffness matrices are $\mathbf{K}(\mathbf{q})_{u}=\mathbf{C}_{u}^{\text{T}}\mathbf{Q}\mathbf{C}_{u}$ and $\mathbf{K}(\mathbf{q})_{s}=\mathbf{C}_{u}^{\text{T}}\mathbf{Q}\mathbf{C}_{s}$; and $\mathbf{Q}=\text{diag}(\mathbf{q})$.
The submatrices $\mathbf{C}_{u}$ and $\mathbf{C}_{s}$ contain the columns of a connectivity matrix $\mathbf{C}\in\{-1,0,1\}^{M\times N}$ corresponding to the free and the fixed nodes of the bar system, respectively.
The connectivity matrix establishes the relationship between connected bars and nodes.
The loads applied to the nodes $\mathbf{P}\in\mathbb{R}^{N\times 3}$ and the positions of the node anchors $\mathbf{X}_{s}\in\mathbb{R}^{N_{s}\times 3}$
are the boundary conditions $\mathbf{b}$ of the form-finding process, and are inputs of the simulator.
To create the boundary conditions vector $\mathbf{b}$, we concatenate $\mathbf{P}$ and $\mathbf{X}_{s}$, and then flatten the resulting matrix into one vector.
Once we know $\mathbf{X}_{u}$, we concatenate it with $\mathbf{X}_{s}$ to obtain all the nodal positions in equilibrium $\mathbf{X}$ of the bar structure.

## Appendix E Implementation details

We implement our work in JAX [Bradbury et al., ] and Equinox [Kidger and Garcia, ].
We use JAX FDM [Pastrana et al., ] as our differentiable simulator and Optax [DeepMind et al., ] for derivatives processing.
For our model and the neural baselines (NN and PINN), we train a separate model instance per task to minimize the loss expectation over a batch of samples of size $B$.
For optimization, we directly minimize the loss for every shape in the batch, one shape at a time.

Training and inference for all models is executed on a CPU, on a Macbook Pro laptop with an M2 chip.
This choice of hardware is motivated by the type of devices structural designers have access to in practice, but also demonstrates the potential for even better model performance on other more powerful hardware, like GPU accelerators.
We train all three models using Adam with default parameters, except for the learning rate, for a fixed number of 10,000 steps.
All the models’ hyperparameters are calibrated to maximize predictive performance on the test set via a random search over 3 different seeds.

We employ SLSQP [Kraft, ] and L-BFGS-B [Zhu et al., ], two deterministic gradient-based optimizers, for direct optimization, as implemented in JAXOPT [Blondel et al., ] with default settings.
To make a fair comparison, we also perform all the direct optimization experiments on CPU and on the same laptop.
Additionally, we impose box constraints on the bar stiffnesses $\mathbf{q}$ (the optimization variables) to prescribe bounds on the stress signs and values as we do for the last layer of our encoder.
We run the optimizers for at most 5,000 steps and with a convergence tolerance of $1\times 10^{-6}$.

Table 3 summarizes the configuration of the neural models studied herein.
We utilize MLPs for the encoder and decoder components of the three models.
The MLPs are a sequence of linear layers followed by a nonlinear activation function.
The weights and biases per layer are initialized at random from a uniform distribution in the interval $[-1/\sqrt{\rho},1/\sqrt{\rho})$, where $\rho$ is equal to the layer’s input size.
We use an ELU [Clevert et al., ] as the activation function in every MLP layer except for the last layer of each encoder.
The encoders employ Softplus as the last nonlinearity to satisfy the strictly positive requirement of our mechanical simulator established in Section 2.3.

The MLP decoder in the NN and PINN baselines takes as input the bar stiffnesses predicted $\mathbf{q}\in\mathbb{R}^{M}$ by the encoder, in addition to a vector with the boundary conditions $\mathbf{b}\in\mathbb{R}^{L}$ input to the mechanical simulator so that this learned decoder sees as much of the information the simulator has access to.
This is different from, for example, conventional PINN approaches where the boundary conditions should be learned in addition to the problem physics [Haghighat et al., , Raissi et al., ]; but akin to the approach proposed in other works [Lu et al., , Bastek and Kochmann, , Mai et al., ] where boundary conditions are treated as either hard constraints or imposed as network inputs.
Here, we simplify the decoder’s goal by imposing the boundary conditions by construction.
This is why the boundaries of the shape targets are matched exactly in all of our examples.
For simplicity, the vector $\mathbf{b}$ only comprises the Cartesian coordinates of the $N_{s}$ fixed nodes and the vertical component of the loads applied to the $N$ nodes (i.e., the horizontal components are zero) per target shape.
The vector $\mathbf{b}$ has size $L=3N_{s}+N$.

Similarly, the last linear layer in the MLP decoders outputs a vector of size $3N_{u}$, representing the 3D coordinates of the free nodes of a bar system $\mathbf{X}_{u}$.
This is motivated by the fact that the mechanical simulator solves the form-finding problem only on the free nodes of a bar system with Equation 7, and it then concatenates the known position of the support nodes $\mathbf{X}_{s}$ to assemble the shape matrix $\mathbf{X}$.
Therefore, we replicate this operation with the MLP decoder so that it predicts as much information as the simulator does.
We then concatenate the output of the last layer of the neural decoder with the known position of the fixed anchor nodes.
We describe task-specific implementation details next.

**Table 3: Neural architectures. We show the composition of the two types of autoencoders we use in this work: a fully neural model employed for the NN and PINN models, and ours with the differentiable simulator as an analytical decoder. In this table, Linear(in_size, out_size) denotes a linear layer with learnable weights and biases followed by a nonlinearity; whereas Concat(in_size, out_size) denotes the concatenation of two vectors. Our model is more compact than the neural baselines because it has at least half of the number of trainable parameters. However, our model takes longer to train due to the computational cost of our mechanical simulator.**
| Task | Model | Encoder | Decoder | $\#$ Parameters $\downarrow$ | Train time (s) $\downarrow$ |
| --- | --- | --- | --- | --- | --- |
| Vaults | NN, PINN | [Linear$(3N,H)$; | [Linear$(M+L,H)$; | 535,412 | 85 |
| $2\times$ Linear$(H,H)$; | $2\times$ Linear$(H,H)$; |  |  |  |  |
| Linear$(H,M)$] | Linear$(H,3N_{u})$; |  |  |  |  |
|  | Concat$(3N_{u}+3N_{s},3N)$] |  |  |  |  |
| Ours | [Linear$(3N,H)$; | Simulator$(M+L,3N)$ | 254,900 | 140 |  |
| $2\times$Linear$(H,H)$; |  |  |  |  |  |
| Linear$(H,M)$] |  |  |  |  |  |
| Towers | NN, PINN | [Linear$(3R,H)$; | [Linear$(M+L,H)$; | 1,245,216 | 136 |
| $4\times$ Linear$(H,H)$; | $4\times$ Linear$(H,H)$; |  |  |  |  |
| Linear$(H,M)$] | Linear$(H,3N_{u})$; |  |  |  |  |
|  | Concat$(3N_{u}+3N_{s},3N)$] |  |  |  |  |
| Ours | [Linear$(3R,H)$; | Simulator$(M+L,3N)$ | 468,880 | 810 |  |
| $4\times$Linear$(H,H)$; |  |  |  |  |  |
| Linear$(H,M)$] |  |  |  |  |  |

### E.1 Masonry shells task

For the masonry vault amortization task, we train MLPs with 2 hidden layers, with $H=256$ units each.
The size of the input and output layers varies as per Table 3.
The batch size is $B=64$.
Our model takes 2 minutes and 20 seconds to train.
The fully neural alternatives, in contrast, take 1 minute and 25 seconds.
We utilize SLSQP [Kraft, ] as our baseline for direct optimization with uniformly sampled, random values of the initial bar stiffnesses $\mathbf{q}$.
We train the PINN baseline to minimize a weighted combination of two terms: the shape loss $\mathcal{L}_{\text{shape}}$ (Equation 2) and the physics loss $\mathcal{L}_{\text{physics}}$ (Equation 5):

$$ $\mathcal{L}=\mathcal{L}_{\text{shape}}+\kappa\mathcal{L}_{\text{physics}}$ (8) $$

Different values of $\kappa$ impact the model performance on this bipartite task (i.e., simultaneously matching target shapes and reducing the residual forces $\mathbf{r}_{i}$).
To strengthen the PINN baseline, we tune the value $\kappa$ in five consecutive increments, from $\kappa=10^{-2}$ to $\kappa=10^{2}$, and train a separate instance of the model for each.
We then pick the PINN that achieves the lowest unweighted loss on the test set (i.e., we set $\kappa=10^{0}$ during inference to evaluate performance).
Table 4 shows our results.
We find that employing $\kappa=10^{1}$ during training produces the best-performing PINN baseline, and we use this PINN for comparison with our model and the baselines in Section 4.

**Table 4: PINN model evaluation on the masonry shells task for five different coefficients $\kappa$ applied to the physics loss $\mathcal{L}_{\text{physics}}$ during training. The table reports the unweighted loss values generated by the trained PINN model predictions at inference time. $\kappa=10^{1}$ produces best results.**
| $\kappa$ | $10^{-2}$ | $10^{-1}$ | $10^{0}$ | $10^{1}$ | $10^{2}$ |
| --- | --- | --- | --- | --- | --- |
| $\mathcal{L}_{\text{shape}}\,\downarrow$ | $1.9~{}\pm~{}0.3$ | $1.5~{}\pm~{}0.3$ | $1.5~{}\pm~{}0.5$ | $\mathbf{3.1~{}\pm~{}1.4}$ | $88.3~{}\pm~{}39.9$ |
| $\mathcal{L}_{\text{physics}}\,\downarrow$ | $33.7~{}\pm~{}12.7$ | $6.2~{}\pm~{}2.5$ | $2.4~{}\pm~{}1.1$ | $\mathbf{0.7~{}\pm~{}0.3}$ | $0.3~{}\pm~{}0.2$ |

### E.2 Cable-net towers task

In the cable-net task, we train MLPs with 4 hidden layers of size $H=256$.
In all models, we reduce the size of the encoder input from $3N$ to $3R$, where $N$ is the total number of nodes in the structure, and $R=48$ corresponds to the number of nodes on the three rings (bottom, middle, and top) that parametrize this task.
This choice reduces the total number of parameters in the encoders and makes the encoder less computationally intensive.
We set the batch size to $B=16$ for all models.
The shape approximation task is underspecified because it only prescribes target height values for the intermediary tensile rings in a given cable-net tower instead of target positions. Therefore, we mask (i.e., multiply by zero) the predicted $x$ and $y$ coordinates of the nodes on these intermediary rings before evaluating the shape loss $\mathcal{L}_{\text{shape}}$.

As for our model, we train it in two stages, reducing the optimizer’s learning rate from one stage to the next.
In the first stage, we optimize for 5,000 steps with a constant learning rate of $1\times~{}10^{-3}$.
In the second stage, we fine-tune our model for another 5,000 steps with a smaller learning rate of $1\times~{}10^{-4}$.
During both stages, the global gradient clip value is 0.01.
The total training time of our model is 13.5 minutes.
We train the MLPs of the NN and the PINN with a constant learning rate of $1\times~{}10^{-3}$ over 10,000 steps.
The training time of both models is equal to 2 minutes and 16 seconds.
The loss function we utilize to train the PINN baseline includes an explicit regularization term, in addition to the shape and the physics losses:

$$ $\mathcal{L}=\mathcal{L}_{\text{shape}}+\kappa\mathcal{L}_{\text{physics}}+ \lambda\mathcal{L}_{\text{reg}}$ (9) $$

Like in the masonry task, we tune the value of the weight coefficient $\kappa$ during training in five distinct steps between $\kappa=10^{-2}$ and $\kappa=10^{2}$.
The value of $\lambda$ is reported in Section 4.
As shown in Table 5, training the PINN with $\kappa=10^{0}$ yields the best performance during training and inference.
We use this PINN variant for comparison with our model and the other baselines.
Lastly, we employ L-BFGS-B [Zhu et al., ], a quasi-Newton method, as our baseline for direct optimization.
We initialize the bar stiffnesses in this task with three different configurations (randomized, expert, and with our trained model) as described in Table 2.

**Table 5: PINN model evaluation on the cable-net tower task for different coefficients $\kappa$ on the physics loss $\mathcal{L}_{\text{physics}}$. The table reports unweighted loss values generated by the trained PINN predictions at inference time. Using $\kappa=10^{0}$ during training produces the best performing PINN model.**
| $\kappa$ | $10^{-2}$ | $10^{-1}$ | $10^{0}$ | $10^{1}$ | $10^{2}$ |
| --- | --- | --- | --- | --- | --- |
| $\mathcal{L}_{\text{shape}}\,\downarrow$ | $42.8~{}\pm~{}17.8$ | $9.6~{}\pm~{}5.6$ | $\mathbf{7.4~{}\pm~{}3.4}$ | $4.2~{}\pm~{}2.6$ | $5.1~{}\pm~{}3.5$ |
| $\mathcal{L}_{\text{physics}}\,\downarrow$ | $1.1~{}\pm~{}0.1$ | $2.1~{}\pm~{}0.0$ | $\mathbf{1.2~{}\pm~{}0.1}$ | $3.8~{}\pm~{}0.7$ | $4.3~{}\pm~{}0.6$ |

## Appendix F Data generation on Bezier surfaces

A Bezier patch $\mathcal{B}$ maps a matrix $\mathbf{C}\in\mathbb{R}^{C\times 3}$ of control points to a smooth surface in $\mathbb{R}^{3}$, $\mathcal{B}:~{}\mathbb{R}^{C\times 3}\to\mathcal{S}(u,v)$, parameterized by local coordinates $(u,v)$.
This parametrization offers clear control over architectural design intent as it enables the exploration of a wide array of smooth geometries by simply changing the positions of a coarse control grid.

### F.1 Data generation for training

A summary of the data generation process is given in Figure 24.
To generate a family of shapes for the shells task in Section 4, we focus on doubly-symmetric shapes generated by a square grid $w=10$ units wide centered on the origin. The grid contains $C=16$ control points arranged in a $4\times 4$ layout.

We then follow three main steps.
First, we vary the 3D coordinates of control points $\mathbf{c}_{1}$ to $\mathbf{c}_{3}$ on a quarter of the control grid.
This construction assures the double symmetry in the generated data.
To vary the position of each of the three control points $\mathbf{c}$, we first sample a translation vector $\mathbf{t}$ at random from a uniform distribution in the interval $[\mathbf{t}_{\text{min}},\mathbf{t}_{\text{max}})$; and we add it to the control point’s reference coordinates $\mathbf{c}_{0}$, such that $\mathbf{c}=\mathbf{c}_{0}+\mathbf{t}$.
The position of $\mathbf{c}_{4}$ is static.
We detail the reference coordinates and the intervals for each control point in Table 6.

Next, we mirror these $4$ control points on the YZ and the XZ Cartesian planes to obtain the position of the remaining $12$ control points in the grid (see callout in Figure 24).
The bounding box of the resulting design space illustrates the scale of the masonry shell task and has dimensions $[2w,2w,w]$.
In the third and last step, we evaluate $N=100$ equally spaced, local coordinates $u$ and $v$ on the Bezier, ranging from 0 to 1, to generate an equal number of points on the surface.
The evaluated points represent the position $\hat{\mathbf{X}}\in\mathbb{R}^{N\times 3}$ of the vertices of the structure we employ as targets to train our model and the baselines.
The coordinates of a point $\hat{\mathbf{x}}\in\mathbb{R}^{3}$ on the Bezier surface are a function of a $(u,v)$ coordinates pair:

$$ $\hat{\mathbf{x}}(u,v)=\sum_{e=1}^{E}\,\sum_{g=1}^{G}\,\gamma_{e}(u)\,\gamma_{g }(v)\,\mathbf{c}_{eg}$ (10) $$

where $\gamma$ denotes a Bernstein polynomial of degree 3; $E=4$ and $G=4$; and $\mathbf{c}_{eg}$ indicates the position of the Bezier’s control points, indexed on a $E\times G$ grid.

Figure: Figure 24: The data generation process of doubly-symmetric shapes for the masonry shells task consists of, first, creating variations of the positions of the points $\mathbf{c}_{1}$ to $\mathbf{c}_{3}$ in a corner of the control grid of a Bezier surface, and, then, mirroring twice. We create a target point $\hat{\mathbf{x}}_{i}$ by evaluating the Bezier at coordinates $(u,v)$.
Refer to caption: extracted/5899952/figures/bezier_surface_points_generation.png

### F.2 Interpolation of Beziers

**Table 6: Generation parameters to sample random variations of the 3D coordinates of the points on a quarter of the control grid of a Bezier surface.**
| Control point | Reference position, $\mathbf{c}_{0}$ | Lower bound, $\mathbf{t}_{\text{min}}$ | Upper bound, $\mathbf{t}_{\text{max}}$ |
| --- | --- | --- | --- |
| $\mathbf{c}_{1}$ | $w/6,w/6,0$ | $0,0,w/10$ | $0,0,w$ |
| $\mathbf{c}_{2}$ | $w/2,w/6,0$ | $-w/2,0,0$ | $w/2,0,w/2$ |
| $\mathbf{c}_{3}$ | $w/6,w/2,0$ | $0,-w/2,0$ | $0,w/2,0$ |
| $\mathbf{c}_{4}$ | $w/2,w/2,0$ | $-$ | $-$ |

We utilize linear interpolation to blend between doubly-symmetric and asymmetric Bezier surfaces.
Since the targets $\hat{\mathbf{X}}$ are a function of the control points matrix $\mathbf{C}$, we interpolate between the control points matrix $\mathbf{C}_{\text{sym}}$ of a surface with double symmetry and that of an asymmetric surface $\mathbf{C}_{\text{asym}}$ to create one design

$$ $\mathbf{C}_{\text{interp}}=(1-\delta)\mathbf{C}_{\text{sym}}+\delta\mathbf{C}_ {\text{asym}}$ (11) $$

where $\delta$ is the interpolation factor and $\mathbf{C}_{\text{interp}}$ is the interpolated control points matrix.
We vary all the control points in $\mathbf{C}_{\text{asym}}$ sampling random translation vectors from uniform distributions like in the doubly-symmetric case, except for the four control points at the corners of the Bezier grid whose position also remains fixed.
We finally generate target points $\hat{\mathbf{X}}$ from $\mathbf{C}_{\text{interp}}$ with Equation 10.
An example of the shapes resulting from interpolating two Bezier surfaces is provided in Figure 19(a).

## Peer Reviews

confidence: 4
contribution: 2
initial_rating: 6
presentation: 2
questions: Overall, the paper is technically sound and proposes an interesting integration of ML and physics. However, there are some concerns as follows:

1. The proposed hybrid approach requires running the physics model during inference, likely increasing training costs. While the authors report inference time, the training time is not provided. Given the limited representational capacity of the network, it must be retrained whenever the design is re-parameterized (a likely scenario during conceptualization). 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.

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.

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.

4. The case studies considered here have relatively low DOFs, and the physics relies on a linear FDM. The authors could discuss the viability of this approach for structures where linear FDM is inapplicable or for dynamic scenarios. In such cases—and even when FDM is applicable but the structure has much higher DOFs—would this approach remain feasible?

5. Furthermore, the variation in MLP and PINN inference times is puzzling. Since the input sizes remain constant, one would expect the inference times for the fully trained networks to be similar. Could the authors comment on this discrepancy?
soundness: 3
strengths: 1. The idea of combining physics with ML for architectural design is promising, as it removes the need for generating costly training labels. Instead, the model can learn by integrating the physics model with inexpensive loss functions.
2. The out-of-distribution performance is also interesting, as it potentially reduces the need for extensive variability in the training data.
summary: The paper introduces a physics-in-the-loop scheme for the design of architectural structures. The authors here use a neural network(MLP) to learn the mapping from the desired structural shape to intermediate mechanical properties (bar stiffnesses) and then use a physics model to generate a physically feasible shape approximating the target.  The authors then apply this method to the design of masonry shells and cable towers,  comparing it with two neural network baselines: one trained to produce physically feasible shapes and the other trained to produce feasible shapes while also ensuring mechanical stability. Results on both these case studies show the proposed approach outperforming both these baselines and being on par with numerical optimization.
weaknesses: 1. The case studies provided are relatively simplistic and do not reflect real-world applications. The physics model used is also quite basic, limiting the method's applicability in practical scenarios.
2. Additionally, as the authors acknowledge, even when trained, the proposed parameterization lacks flexibility and requires retraining whenever the design representation changes.

---

confidence: 5
contribution: 2
initial_rating: 5
presentation: 3
questions: 1. Does the term *fabrication requirements* imply that each bar in the shell examples has to follow additional constraints?
2. Provide more discussion on the usage of $\tau$ as pointed out in the previous section
3. Revise the manuscript to discuss the decoder in the main body of the manuscript.

Minor issues:
1. Define (or discuss) *mechanical integrity* and *mechanical efficiency* earlier in the manuscript. What do they exactly mean?
soundness: 3
strengths: The paper is nicely written, although there are a few small details which could make the paper more understandable for neurips attendees, who may not have prior knowledge on the mechanical structure presented.

The method presented appears sound. The architecture as shown in Fig 1, leveraging JAX, makes intuitive sense.
summary: This manuscript presents an interesting application of AE-type of network architecture to facilitate design and optimization of mechanical structures, mostly in (building) architectural structures. Although the method presented is mostly ad-hoc, the manuscript is well written, the proposed method appears sound, and the potential application is clearly indicated.
weaknesses: Concerns to the paper:

1. Although it is not clearly explained, it appears that each trained model is closely related to the underlying architectural structure (masonry shells or cable-net tower). The question is whether the trained model is generalizable. Even for the shells, if $M$ (number of bars) and $N$ (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. 

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.

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. 

4. Structure of the paper: the decoder is not discussed in the main body of the paper, mostly in the appendices. This need to be fixed. 

Minor issues:
The terms *mechanical integrity* and *mechanical efficiency*  should be explained early in the paper to make it easier for readers to appreciate the paper.

---

confidence: 2
contribution: 2
initial_rating: 5
presentation: 3
questions: * Are the authors going to be releasing the code? As an application paper, it would be beneficial to the community if they were.
* Does the chosen bar stiffness affect the loading on the design since a stiffer bar might lead to a heavier design? 
* Could you provide further clarification on how the stress direction is determined a-priori? 
* 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.
* What mechanism is in place to prevent the model providing excessively large stiffness values for all the bars?
* 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)?
soundness: 3
strengths: * The presentation of the approach is nicely structured, and the figures are clean and easy to understand. 
* As an application paper, the work demonstrates a nice implementation of combining a differentiable simulator (as the decoder) with a neural network (as the encoder). The authors also build a physical prototype (Figure 6) based on the output of the approach. It is highly unusual to take this extra step to build a physical prototype and should therefore be commended.
* The results successfully highlight how using the simulator ensures that the final design results in a pin-jointed structure where the resolved forces at the internal pins are at equilibrium.
summary: This paper combines neural networks with a differentiable mechanics simulator to develop a model that accelerates the solution of shape approximation problems for architectural structures modeled as bar systems. The proposed approach achieves better accuracy and generalization than the neural network approach and gets comparable accuracy to direct optimization. The two design domains are masonry shells and cable-net towers.
weaknesses: * The motivation for this approach is not entirely clear. The design pipeline appears to start with a $\tilde{\mathbf{X}}$ and then produces a set of stiffness values. These stiffness values, along with the boundary conditions are then fed to the mechanical simulator to result in the final design (presumably, with the neural network’s stiffness values):
    * In practice, how does one come up with a design $\tilde{\mathbf{X}}$? 
    * How does the simulator only require $\mathbf{p}$ and $\mathbf{b}$ for shape prediction? It seems like the simulator is missing information to build a shell with $N$ pins. 
* The cost of amortization (i.e. the number of direct optimization runs to build the data set for training) does not seem to be mentioned as a limitation of the approach. Are there instances in the real world where amortization of the optimization routine would be needed? For example, with a training set, it seems like 100 x direct optimization is required to train the model. This overhead would need at least 100 runs of the proposed approach to make it well-motivated. In a real design pipeline, how often would it be expected to call such a function? It already runs in less than 2 seconds for cable-net towers and around 6 seconds for the masonry shell. Is this this wall-clock time a problem in practice? If so, the paper would benefit from this written explicitly. 

Minor
* $p$ for norm and $\mathbf{p}$ for loads are a bit confusing in terms of notation.
* Suggestion: use test-time inference rather than inference for predictions, since inference can also refer to inference over parameters.

---

confidence: 3
contribution: 3
initial_rating: 8
presentation: 3
questions: See weakness section
soundness: 3
strengths: The paper presentation is very good, the motivation is clear, and the method is justified well. Experimental results are extensive. The authors also provide a real life sample designed from their algorithm, which the readers will find very interesting.
summary: The authors propose an AI-driven design algorithm for architectural structure that can generate mechanically feasible designs for bar systems according to a target geometry. Previous approaches such as pure neural networks cannot guarantee mechanical design requirements. Traditional direct optimization methods can satisfy mechanical criteria, however are inefficient. The authors propose a hybrid scheme, combining neural networks and mechanical simulators to achieve computational efficiency and satisfaction of mechanical constraints. The authors approach consists of a encoder-decoder neural network, where the decoder is replaced by a differentiable mechanical simulator. The encoder structure provides efficiency, while the decoder ensures the design satisfies mechanical constraints. Experimental results show that the authors method achieves satisfactory performance in terms of design and computational efficiency.
weaknesses: The authors could do some revision of the figure. In figure 1, it seems that 3 figures are simultaneously provided to the encoder network, although the text description indicates only one shape is provided. Please take care to avoid such confusion from readers. Also describe the dimension of the matrix K. In figure 3, it is hard to see how the x vectors are being updated. Zooming in and adding some texts will be better for understanding.

Meta-review: The paper introduces a surrogate model for designing architectural shape structures modeled as bar systems, a simplified structural representation consisting of nodes connected by bars. This model combines a neural network (NN) with a differentiable mechanical simulator. Given an input target shape, provided, for example, by a designer, the NN computes stiffness parameters that encode the mechanical behavior of the structure, while the simulator outputs the corresponding structure. The model is trained end-to-end. Validation is performed on two problems: masonry shells and cable-net towers. The baselines consist of an encoder-decoder NN, the same architecture augmented with an additional mechanical loss (similar to a physics-informed neural networks (PINNs) approach), and a traditional direct optimization. The proposed hybrid model approaches the performance of direct simulation methods at a much lower computational cost.

This is an application paper, and the idea of combining differentiable physical models with NNs has been explored in various fields. The originality of this work lies in its application domain. The reviewers consider this work, even though it is evaluated on relatively simple mechanical design problems, to be a significant first step toward developing practical applications. The authors’ responses and the additional experiments provided during the rebuttal adequately address the reviewers' concerns. I recommend accepting the paper.

Decision: Accept (Poster)