## 3 Computation of Minimum Thickness  
The numerical formulation to compute the spatial configuration of thrust networks is based on the equilibrium with independent force densities proposed in [21]. The horizontal projection of the considered networks is kept fixed and this projection is called _form diagram_[20]. With such a formulation, the spatial equilibrium of the networks can be computed in terms of two sets of parameters: the independent force densities (\(\mathbf{q}_{\text{ind}}\)) and the heights of the supports (\(\mathbf{z}_{\text{b}}\)). From the vertical equilibrium, the height of any internal node of the network \(\mathbf{z}_{\text{i}}\) can be calculated as:  
\[\mathbf{z}_{\text{i}}=(\mathbf{C}_{\text{i}}^{\text{T}}\mathbf{Q}\mathbf{C}_{ \text{i}})^{-1}(\mathbf{p}_{\text{z}}-(\mathbf{C}_{\text{i}}^{\text{T}} \mathbf{Q}\mathbf{C}_{\text{b}})\mathbf{z}_{\text{b}}), \tag{2}\]  
where \(\mathbf{C}\) represents the connectivity matrix of the network, defined after [38]. \(\mathbf{C}_{\text{i}}\) and \(\mathbf{C}_{\text{b}}\) are slices of the connectivity matrix corresponding to the internal and boundary/restrained nodes, respectively. The vector \(\mathbf{p}_{\text{z}}\) collects the self weight lumped to the internal nodes of the network (i.e., using the tributary area of each vertex with respect to the middle surface, as illustrated in Section 4). \(\mathbf{Q}=\text{diag}(\mathbf{q})\), where the vector of force densities \(\mathbf{q}\) accounts for both dependent (\(\mathbf{q}_{\text{dep}}\)) and independent (\(\mathbf{q}_{\text{ind}}\)) force densities which are linked by the following horizontal equilibrium equation:  
\[\mathbf{q}_{\text{dep}}=-\mathbf{E}_{\text{dep}}^{-1}\,\mathbf{E}_{\text{ind }}\mathbf{q}_{\text{ind}}, \tag{3}\]where \(\mathbf{E}\) is the equilibrium matrix computed as in [21], and \(\mathbf{E}_{\text{dep}}\) and \(\mathbf{E}_{\text{ind}}\) are slices of the equilibrium matrix considering the dependent and independent force densities, respectively. This mathematical formulation is general; i.e., it is not limited to any specific geometry or support condition, and could be used to model different three-dimensional systems, modelled with different form diagrams. More details about this numerical formulation and the selection of independent edges can be found in [25]. From this numerical formulation, the minimum thickness is calculated as a constrained nonlinear optimisation problem (NLP) considering the analytical description of the vault's upper and lower bounds, i.e., extrados and intrados. The optimisation problem is described below:  
\[\underset{\mathbf{q}_{\text{ind}},\mathbf{z}_{\text{b}},t}{ \text{minimise}} t,\] (4a) subject to \[\mathbf{q}_{\text{dep}}\geq\mathbf{0}, \tag{4b}\] \[z_{i}^{\text{LB}}(t)\leq z_{i}(\mathbf{q}_{\text{ind}},\mathbf{ z}_{\text{b}})\leq z_{i}^{\text{UB}}(t),\] (4c) \[t\geq 0. \tag{4d}\]  
The linear objective function (4a) considers the structural thickness \(t\) only. Constraints (4b) represent the compression-only requirement from [PERSON]'s assumptions. Constraints (4c) enforce the network to be contained within the bounds of the structure, calculated for the thickness \(t\). Constraint (4d) expresses that the thickness must be positive. The gradients and Jacobian matrix of the problem are computed analytically. For an exhaustive discussion about this NLP formulation, the reader is referred to [25]. To solve this NLP, a combination of two solvers was applied, combining a sequential least squares programming (SLSQP) solver [39], and an interior-point optimisation (IPOPT) solver [40]. Both are suitable to solve constrained NLP, the latter being a more robust solver, which is used to polish the results, and to deal with the problems with larger numbers of variables (e.g., more than 100 variables). The objective function (4a) can be accordingly modified in order to minimise, or maximise the (horizontal) thrust of the structure, once the thickness has been fixed.