From Stress to Shape: Equilibrium of Cloister and Cross Vaults

: The assessment of the equilibrium and the safety of masonry vaults is of high relevance for the conservation and restoration of historical heritage. In the literature many approaches have been proposed for this tasks, starting from the 17th century. In this work we focus on the Membrane Equilibrium Analysis, developed under the Heyman’s theory of Limit Analysis. Within this theory, the equilibrium of a vault is assessed if it is possible to ﬁnd at least one membrane surface, between the volume of the vaults, being in equilibrium under the given loads through a purely compressive stress ﬁeld. The equilibrium of membranes is described by a second order partial differential equation, which is deﬁnitely elliptic only when a negative semideﬁnite stress is assigned, and the shape is the unknown of the problem. The proposed algorithm aims at ﬁnding membrane shapes, entirely comprised between the geometry of the vault, in equilibrium with admissible stress ﬁelds, through the minimization of an error function with respect to shape parameters of the stress potential, and then, with respect to the boundary values of the membrane shape. The application to two test cases shows the viability of this tool for the assessment of the equilibrium of existing masonry vaults.


Introduction
The assessment of the equilibrium of masonry structure is a topic of high relevance, given the diffuse presence of masonry buildings in Italy and Europe, often with high historical value.
Over the centuries, different methodologies have been developed for the assessment of the equilibrium of masonry constructions and to quantify their degree of safety, starting from practical rules [1], until the definition of graphical methods applied to two dimensional structures (such as portals, bridges, and flying buttresses [2,3]) and to axis-symmetric vaults under monoaxial stress regimes, such as masonry domes [4]. Graphical methods have been extended to consider also biaxial stress states, again under the hypothesis of axial symmetry of the structure [5].
For masonry material, the starting hypotheses have been clearly formulated by Heyman in [6,7], and assume that the material is infinitely resistant in compression and that sliding between rigid blocks does not occur inside the material. Under these assumptions, the theorems of Limit Analysis hold, and an equilibrium solution is defined as admissible if it is balanced with the given loads and is everywhere of pure compression. In particular, based on the Safe Theorem, the impossibility of collapse of a masonry structure, under the given loads, is ensured if at least one admissible equilibrium solution can be found.
A wide number of methods have been developed for the assessment of masonry vaults within the Heyman's theory, both considering the kinematic and the Static Theorem. Concerning the Kinematic theorem of Limit Analysis, the proposed methods are based on energy criteria minimization. Indeed, the minimum of the total potential energy can be used to look at both stable and non-stable mechanisms [8,9]. Moreover, it can be also used to couple internal stress states with the corresponding mechanisms as in [10]. Other methods, such as the Thrust Network Analysis (TNA), are based on the application of the Static Theorem. TNA, introduced in [11] and developed in [12][13][14][15][16], consists in the discretization of a masonry vault through a finite number of compressed bars, in equilibrium with the load of the vault, and located in-between the intrados and the extrados [17].
Among different continuum approach methods, here we focus on the Membrane Equilibrium Analysis (MEA), first proposed in [18], and further developed in [19]. For the assessment of the equilibrium of a vault, MEA requires to find at least one membrane surface comprised in-between the geometrical boundary of the vault, balancing the external loads through a negative semi-definite membrane stress field.
In MEA, the equilibrium is reduced to a scalar second order partial differential equation with non-constant coefficients. If the shape of the membrane is assigned, the coefficients are the components of the curvature of the membrane, and the unknown of the problem is the stress potential of the membrane stress projected onto a reference plane (Pucher stress). On the contrary, if the stress potential is assigned, the coefficients of the equation are the components of the Hessian of the stress potential and the unknown is the scalar function defining the elevation of the surface with respect to the reference plane. In both cases, the equilibrium equation may change its behaviour (hyperbolic-parabolicelliptic) depending on the assigned shape or on the assigned stress potential. This has a strong implication on the numerical solution of the equilibrium equation, that can be solved as a boundary value problem when the equation is elliptic, or as an evolutionary problem when the equilibrium equation is hyperbolic at least in a part of the domain.
MEA has been adopted mainly for the assessment of the equilibrium of masonry constructions [20][21][22][23][24], but it is a general tool for computing the stress field on curved membrane surfaces. For example, it has been used for the evaluation of the stress field of soft tissues, such as the human cornea [25].
An approach similar to MEA, based on the search of the Airy stress function, has been proposed in [26], where NURBS are used for the discretization, exploiting their high-order continuity.
In [27] a numerical iterative procedure for the assessment of the equilibrium of masonry vaults, based on LSM (a FEM-like approximation of the governing equation, see [28]), has been introduced. Starting from a tentative membrane shape, included in the volume of the vault, the method looks for an equilibrated concave stress potential, ensuring the negative semi-definiteness of the corresponding stress. This procedure, however, is not easily applied to non-concave shapes, such as cross vaults given the non-ellipticity of the corresponding equilibrium equation. Recently, in [29] a numerical procedure based on standard FEM, for the design of purely compressed shells, has been put forward. With this method, the starting point is a class of concave stress potentials, and the equilibrium equation is purely elliptic.
In the present paper, focusing on the assessment of the equilibrium of existing masonry vaults, we start from a family of concave stress potentials in order to ensure the admissibility of the stress field. This implies the ellipticity of the equilibrium equation: Therefore, we look for an optimal membrane shape comprised between the geometrical boundary of the vault, as required by the Safe Theorem of Limit Analysis. We approach this problem through an optimization algorithm, minimizing an error function first with respect to the parameters of the stress potential, and then with respect to the boundary conditions of the the membrane shape.
Optimization is a recurring tool in different engineering applications. In the present work, we find an optimal ideal membrane shape, which best satisfies the requirements of the Safe Theorem of Limit Analysis, which result in some mathematical constraints. We point out that this is not the case of a shape optimization, since the shape of the analyzed vaults is fixed, and we just optimize an ideal membrane. In other engineering application, however, the shape of structures is optimized under the constraints of material saving, or compliance minimization. This is the case of the Structural Topology Optimization [30], used for optimally distribute the material to minimize the compliance o structures under given loads [31][32][33]. Optimization is also used in geotechnical applications [34], optimization of geodesic domes [35], of soil-steel structures [36], and lightweight structure [37]

Mathematical Preliminaries
In this section some mathematical details concerning the formulation of the equilibrium equation of a membrane are recalled, the notation is introduced, and the unilateral constraints on the material are introduced.

Geometry and Equilibrium
Referring to Figure 1 for notation, the geometry of a generic membrane surface S contained inside a typical shell structure can be described in Monge form in terms of the position vector x of the points of S as follows: where {e 1 , e 2 , e 3 } is the orthonormal triad coherent with a given Cartesian frame {O; x 1 , x 2 , x 3 }, the couple of coordinates (x 1 , x 2 ), which are the parameters of the curvilinear description (1), belongs to a region Ω of the plane {O; x 1 , x 2 } called planform of S, and f = f (x 1 , x 2 ) is the rise of the membrane with respect to the x 1 − x 2 plane. We assume f ∈ C 0 (Ω)-that is to say that the considered surface S is continuous but not necessarily smooth. An efficient way to describe the membrane equilibrium of a membrane under a surface load q per unit area, is the formulation proposed by Pucher in [38]. Pucher analysis is based on the introduction of the projected stress S, in terms of which two of the three equilibrium equations for the membrane can be made independent of the membrane shape.
The membrane stresses are tangent in each point to the membrane surface S. Therefore it is convenient to introduce a curvilinear basis on the membrane, which are computed as In addition, it is useful to introduce the contravariant basis a α . α = 1, 2 defined such that a α · a β = δ β α , being δ β α the Kronecker delta.
Therefore we can decompose the membrane generalized stress tensor as where T αβ are the contravariant components of the generalized stress. The projected (Pucher) stress components S αβ (see Figure 1) are easily defined in terms of T αβ as S αβ = J T αβ , J being the ratio between the surface area and the projected area. In terms of Pucher stress components, the vector equilibrium equation, projected into the nonorthonormal vector basis e 1 , e 2 , n = a 3 , n being the unit normal to S, becomes: where {p 1 , p 2 , p 3 } are the Cartesian components in the given reference of p = J q, that is the load per unit projected area. By using the projected stress, the first two equilibrium equations do not depend on the membrane shape and are the same as for the plane stress problem. In the case of pure vertical loading, say p = pe 3 , the first two equations can be solved introducing an Airy stress potential F(x 1 , x 2 )-that here we only assume to be continuous-defined such that The third equation of (4) expresses the balance of the vertical component of the force p 3 = p · e 3 with the scalar product between the Pucher stress tensor and the Hessian of the function f , that defines the curvature of the surface. On introducing (5) into (4), we obtiain, as transversal equilibrium equation Equation (6) is a second order partial differential equation, that can be either elliptic, parabolic or hyperbolic, depending on its (usually non-constant) coefficients. If the problem is hyperbolic or parabolic, the characteristic lines are real and the problem must be formulated as an evolutionary problems by assigning two scalar conditions on a non-characteristic line.
If Equation (6) is elliptic, the problem can be formulated as a boundary value problem and the boundary conditions for this equation can be either of the Dirichlet type, i.e., or Neumann type, i.e., ∂F ∂n where p and q are known function on the boundary ∂Ω, and n is the outward unit vector normal to the boundary. The advantage of boundary value problems is that the solution of the corresponding discretized partial differential equation can be obtained by a simple matrix inversion.

Unilateral Constrains and Singular Stresses
The restrictions that the generalized stress T is negative semi-definite and does not work for the corresponding strain E, that is positive semi-definite, define Rigid No-Tension (RNT) materials in the sense of Heyman. In formulae: The first application of Pucher's transformation for NT masonry vaults can be found in [18], where it is shown that, due to the NT constraint, both the surface generalized stress and the matrix of the projected stresses must be negative semi-definite. In terms of the stress function F, this condition can be written as The condition of semi-definiteness of the stress correspond to concavity of the stress potential.

Outline of the Method
The Membrane Equilibrium Analysis (MEA) applied to masonry vaults, requires to find a membrane surface located in between the intrados and the extrados of the vault and, according to Heyman's assumptions, sustaining the loads through compressive stresses. The membrane is a determined structure, and this last constraint is generally not satisfied on a given membrane shape for given loads. In such a case, the membrane shape has to be changed iteratively from an initial guess shape, until one that satisfies the unilateral restrictions on the stress is found.

Unknowns and Data in the Governing Equilibrium Equation
On adopting a numerical discretization method based on a variational approximation of the governing Equation (6), such as the Finite Element Method, or the Lumped Stress Method defined below, the differential equation must be elliptic. Therefore, if the unknown is the stress and the given form of the membrane is not everywhere convex or concave, these methods of approximation fail and cannot be applied. A possible way out, already proposed in [29], consists in switching the roles of unknowns and data, by assuming a parametric family of concave stress potentials as given, and considering the shape as the unknown. Then, the equation is elliptic for any choice of the parameters defining a specific member of the given family of stress potentials, and the corresponding shape can be obtained by prescribing convenient boundary data for the shape. Besides, for any such concave stress potential, there exists a small repertoire of purely compressed equilibrium shapes, that one can also explore, depending on the possibly different boundary data, that can be imposed on the shape. In the present paper, the equilibrium solution is tackled by constructing an LSM approximation of the fundamental second order differential equation of equilibrium (6). This computer routine, which returns a shape associated to any given stress function, can be manipulated by changing, following a minimization algorithm, the parameters controlling the stress potential and the boundary data for the shape.

Numerical Discretization and Algorithm
The transverse equilibrium Equation (6) is solved numerically by discretizing the domain into a finite number of triangular elements, and adopting a Lumped Stress Method [28] approximation of the Hessian of both the stress potential F and of the membrane shape f . With the LSM approximation, the two functions are described by simplicial surfaces, defined by the nodal valuesF i andf i , with i = 1, 2, . . . , N, N being the total number of nodes of the given triangulation.
The vectors collecting the nodal values of the shape function f and of the stress potential F are denotedf andF.f is divided into two vectorsf I andf B , of length N I and N B , distinguishing between internal and boundary nodes, respectively.
The search of a statically admissible stress regime for the vault is carried out through an algorithm divided in several steps. Given the geometry of the vault, defined by its intrados and extrados, we start on assigning an appropriate functional form for the stress potential, depending on a limited number n of parameters, such that F = F(x 1 , x 2 ; a 1 , a 2 , . . . , a n ).
A first step is the parametric solution of the discretized version of Equation (6), with suitable boundary conditions, in terms of the membrane shape defined by f ; therefore, f = f (x 1 , x 2 ; a 1 , a 2 , . . . , a n ).
The second step consists in the optimization of the parameters a 1 , a 2 , . . . , a n , performed by minimizing an error function I defined as the integral of the quadratic approximation of the characteristic function of the vault domain.
Given a point x ∈ Ω, the characteristic function indicates whether the corresponding values of the rise of the membrane is comprised between the volume of the vault, or not.
In particular Its quadratic approximation reads For a generic point x of the domain, given f intr (x) and f extr (x), the approximated characteristic function is represented in Figure 2 f(x) We denote Ω c the subregion of the entire planform Ω obtained by removing from the domain a boundary strip of size c, defined as a fraction of the diameter L of Ω (in the applications we set c = L/5). The rationale behind this procedure is linked to the nature of the governing equation: due to ellipticity, the solution is only weakly dependent on the boundary data, that is, altering the boundary values has a detectable effect only in a narrow band close to the boundary. So, we can split the optimization into an interior and an exterior part, postponing the optimization of the boundary data.
The error function I is then defined as follows: Ω a and Ω b being the parts of the domain Ω c , for which the membrane surface is below the intrados, or above the extrados, respectively. In the third step, the values of a 1 , a 2 , . . . , a n , found in the previous step, are fixed and the boundary values of the discretized membrane shape are considered as variables. Then, the discretized version of Equation (6) is solved parametrically with respect to these boundary values. In the vector notation introduced above, the transverse equilibrium equation reads (K xx,Hjfj )(K yy,HhFh ) + (K yy,Hjfj )K xx,HhFh ) − 2(K xy,Hjfj )(K xy,HhFh ) +p H = 0 where the summation convention on repeated indices is assumed on the lower-case indices (j, h = 1, 2, . . . , N), whilst the upper-case index H(H = 1, 2, . . . , N I ) is not summed, and K xx , K yy , and K xy are coefficient matrices computed via the LSM scheme, representing the discrete counterparts of the second derivatives ∂ 2 /∂x 2 1 , ∂ 2 /∂x 2 2 , and ∂ 2 /∂x 1 ∂x 2 . Equation (15) can be rewritten in a compact form by factoring out the unknown term f, as where K ij = (K yy,HhFh )K xx,Hj + (K xx,HhFh )K yy,Hj − 2(K xy,HhFh )K xy,Hj .
Expression (16) can be also rewritten distinguishing between boundary conditions and domain equations, in the form where K I I and K IB are submatrices of K corresponding to rows and columns associated to internal and boundary nodes. Finally,f B is the vector of the boundary parameters. From Equation (18),f I can be evaluated aŝ On introducing (19) into a discrete version I 2 of the objective function I, we get We then optimize I 2 (f I ) with respect tof B . Notice that it is possible to compute explicitly the gradient of I 2 ( f I ( f B )) through the chain rule, in the form that is The membrane shape obtained at the end of this minimization process can be compared to the actual geometry of the vault. If the final value of the objective function I 2 is zero, the membrane is completely comprised between the intrados and the extrados of the considered vault, and therefore the solution is statically admissible and proves the stability of the vault under the given loads, according to the Safe Theorem of Limit Analysis. In such a case, the degree of safeness of the structure can be assessed by reducing the slenderness of the vault and computing a geometrical safety factor, as proposed by Heyman [6].

Cloister Vault
Consider the case of a cloister vault under its self-weight, whose geometry is described by the volume enclosed in-between two surfaces, the intrados and the extrados, defined canonically inside the domain Ω = {(x 1 , x 2 ) ∈ [−L, L] × [−L, L]} as follows: where Ω 1 , Ω 2 , Ω 3 , and Ω 4 are the subregions of Ω depicted in Figure 3.  The values adopted in Equations (23) and (24) are L = 2.5 m, h = 2.13 m, and t = 0.23 m. The values have been taken by the geometry of the vault analyzed in [21], and located in the interior of Palazzo Caracciolo di Avellino in Naples. From [21] we also assume p = 8 kN/m 2 as distributed load per unit projected area.
By following the methodology described in Section 3, the transverse equilibrium Equation (6) can be solved for the geometry f (x 1 , x 2 ), by assuming a special parametric form for the stress potential.
For the case at hand, inspired by the analysis described in [21], we identify a central zone Ω * of the domain Ω, in which the projected stress is a plane pressure. Ω * is delimited by a curve Γ defined as r c (θ) = r 0 where r 0 is a parameter, α is set to 0.05 and P to 40. In Figure 5 we show the contour r c obtained for two representative values of r 0 , namely r 0 = 0.5 and r 0 = 1. Inside Ω * , the Pucher stress is assumed as isotropic, that is The corresponding stress potential is therefore easily computed as where p 0 is a parameter defining the intensity of the stress potential (i.e., of the stress) that will be subsequently optimized. In the external zone Ω ext the stress is assumed to be quasi-uniaxial in the radial direction. This means that, on definingr as the direction of the radii emanating from the geometrical center of the domain and directed toward the external boundary, the stress in Ω ext can be written as S = −σ(s)r ⊗r + εσ 0t ⊗t (28) where σ(s) is related to p 0 and denotes the stress intensity along the radius, and εσ 0 is a small stress which has the function to give a small extra-curvature to the potential in the directiont orthogonal tor, in order to preserve the ellipticity of Equation (6) when solved for the shape f (x 1 , x 2 ). The stress potential associated to the stress distribution (28) in the external zone is obtained by imposing a linear growth in the radial direction, enforcing the continuity of the potential F and of its slope across the curve Γ: where s is a parameter indicating the distance from Γ along the radius. In Figures 6, we show the the stress potential for two representative values of r 0 . Adopting these two shapes as the stress potential, and taking as boundary condition for the membrane geometry f (x 1 , x 2 ) the value of the rise of the middle surface between the intrados and the extrados of the vault, Equation (6) can be solved using LSM, as described in Section 3.2. By constructing a structured triangulation based on a regular node distribution over the planform Ω, the membrane shapes corresponding to the stress potentials of Figure 6, are represented in Figure 7. Figure 7 gives an idea on how the choice of the parameter r 0 (and thus, the extension of the central zone characterized by a biaxial stress state) influences the final shape of the membrane geometry. In the second case (r 0 = 1), the membrane geometry resembles the shape of a cloister vault more than in the case of r 0 = 0.5. It is also evident that the height of the two shapes (in the direction of the x 3 -axis) are not comparable with the height of the cloister vault defined from Equations (23) and (24), which is 2.13 m. The minimization of the error function (20) with respect to the parameters p 0 and r 0 is performed as follows: first, r 0 is given a finite number of values in the range [0, L], i.e., r 0,k = 1.1 + k0.01, k = 0, . . . , 20. Then, for each value of r 0,k , we search the p 0,k that minimizes I k = I k (r 0,k , p 0,k ). In Figure 8 we show the considered values of the r 0,k on the horizontal axis, and on the vertical axes we show the corresponding optimal value of I k .
The lowest value of I k is reached when r k = 1.14, which correspond to p 0,k = 3.45 × 10 3 . Therefore the parameters of the stress potential are fixed to r 0 = 1.14 and p 0 = 3.45 × 10 3 . The best-fitting membrane geometry is then represented in Figure 9.
In Figure 10 we also show two sections of the membrane and of the intrados and extrados, the first one taken at x 2 = 0, and the second one taken along the diagonal of the domain, that is for x 2 = −x 1 .  Figure 8. Plot of I k as function of r k . Figure 9. Best fitting membrane shape after the first optimization phase, for the problem of the cloister vault. We notice from Figure 10b that some points of the membrane fall outside the volume of the vault; this is confirmed also by Figure 11, where we show the color map of the function g(x 1 , x 2 ) defined as (30) for which positive values indicate points in which the membrane shape is above the extrados of the vault, negative values indicate when the membrane is below the intrados, and vanishing values indicate nodes that are comprised within the volume of the vault. We then perform the last step of the optimization, with respect to the boundary values of the membrane shape, obtaining a new membrane geometry whose sections for x 2 = 0 and x 2 = −x 1 are shown in Figure 12. In this case, we notice that the points that were falling outside the volume of the vault before this optimization step, are now comprised between the intrados and the extrados. The whole algorithm has been repeated assuming the same loads, but reduced thickness, in order to assess the geometrical safety factor. The higher thickness reduction for which the proposed algorithm still produces admissible solutions, in the sense of Limit Analysis, is 38% of the original thickness, corresponding to a geometrical safety factor of 1.61.

Cross Vault of the Anterior Portico of the Church of St. Pietro in Vineis in Anagni
In this section we consider the case of the anterior portico of the church of San Pietro in Vineis in Anagni. The portico consists of three similar cross vaults: the geometrical description of each vault, referring for notation to Figures 13 and 14, is given by:  The values given to the geometrical parameters introduced in Equations (31) and (32), are reported in Table 1. Table 1. Intrados and extrados data referring to notation in Figure 14. For this kind of vault a concave membrane surface comprised between the intrados and the extrados does not exist; as a consequence, using Equation (6) with a known shape and unknown stress potential would lead to a non-elliptic problem, that cannot be solved using a FEM-like procedure. The inverse problem, consisting in the search of the shape corresponding to a given concave stress potential, is elliptic instead, and the FEM-like procedure proposed in [21] can be applied.
It is a well-known result (see [27]) that the stress potential corresponding to a standard cross vault has the shape of a cloister vault. The case at hand, however, is not a standard cross vault since there is a curvature not only orthogonally to each median of the planform, but also in the median direction. We expect, therefore, a biaxial stress field, with components both along the direction of the medians, and in its orthogonal direction.
We propose therefore the following stress potential F(x 1 , For sake of example, in Figure 15 we present three different stress potentials obtained for β = 0, β = 0.3, and β = 0.7 and α 0 = 1. We assume the same definition of the error introduced in Equation (20) and evaluate it for different values of β, and then minimizing I k with respect to α 0 . The values of β considered in this case are The corresponding plot of the error I k versus β is reported in Figure 16 in a logarithmic scale along the error axis. The minimum occurs at β = 0.4, and corresponds to an error of 3.94 × 10 −2 . The corresponding stress potential is shown in Figure 17. The corresponding shape of the membrane, obtained through the numerical solution of Equation (6), and imposing as boundary values for f (x 1 , x 2 ) the values of the midsurface (between the intrados and the extrados), is shown in Figure 18.
To evaluate the position of the membrane surface after this first optimization phase, we show in Figure 19 the colormap of g(x 1 , x 2 ), defined like in Equation (30), applied to this case study. The presence in Figure 19 of points for which the error function is non zero indicates that the membrane found after the first optimization step is not fully comprised between the intrados and the extrados. Therefore we need a further optimization step to achieve this goal.  The second step of the minimization process is then performed with respect to the boundary values of the membrane shape. The value of the error function I = I(f B ), obtained as a result of the second optimization phase, is zero everywhere, meaning that the membrane fits in-between the bounding surfaces of the vault.
For sake of example, we show in Figure 20 a graph of a section of the intrados, the extrados and the membrane shape along the diagonal line x 2 = − b i a i x 1 , that, after the first minimization phase, exhibits points lying outside the vault geometry. As expected, there are no points of the membrane section falling outside the geometry of the vault.
The influence of changing the boundary valuesf B is made evident by the color map of the difference h(x 1 , x 2 ) defined as and shown in Figure 21, where f is the membrane geometry at the end of the first minimization step, and f opt is the geometry after the final minimization with respect to the boundary conditions.  From Figure 21 we see that the main effects of the modification of the boundary conditions is that the rise change of the membrane, represented by Equation (35), after the second optimization phase, are visible mainly in the boundary part of the domain, as expected from the discussion about elliptic partial differential equations. Moreover, we see that h(x 1 , x 2 ) is negative in the lowest-right and upper-left corners, meaning that the membrane rise is diminished after this optimization, fixing the error shown in Figure 19.
In order to quantify the Geometrical Safety Factor for this vault, the algorithm has been repeated progressively reducing the thickness. The maximum possible reduction, given the assumed shape for the stress potential, is 36%, corresponding to a GSF= 1.56.

Discussion
Concerning the algorithm adopted in both examples, we point out that it is an automatic procedure that is easily adaptable to the specific case. The only different step between the proposed algorithms, beyond the geometries of the vault, is the parametric description of the stress potential, which is specific to each case study, and consists, from an operative point of view, in replacing a computer routine to another one.
Concerning the Geometrical Safety Factor, computed at the end of each application, we remark thaat it is an intrinsic characteristic of each structure, and does not depend, in general, on the adopted membrane shape or its associated equilibrated stress. With computational models, such as MEA, it is only possible to try to find the best value, that, in turn, is associated to specific choices on the stress potential (i.e., its shape, or its parametric description), or on the discretization, or on the specific algorithm.
In our examples, we made two precise choices for the shape of the stress potential, based on experience and on literature examples. A richer description of the stress potential could have improved the obtained value of the Geometrical Safety Factor, at the price of a more complicate and slow optimization algorithm, especially for the first optimization pahse, which is characterized, for each evaluation of the objective function, by a matrix inversion, whose computational cost, in turn, depends in the adopted discretization. Therefore, the proposed values for the GSF have to be considered as conservative values, and could be improved by using more refined optimization algorithms, not yet developed in the framework of the Membrane Equilibrium Analysis.

Conclusions
In this work we present a method for the assessment of the equilibrium of masonry vaults based on the search for a suitable membrane shape in equilibrium with a given compressive stress field. The membrane shape is required to be entirely comprised between the intrados and the extrados of a vault, according to the Safe Theorem of Limit Analysis.
The proposed method aims at finding a concave stress potential in equilibrium with the applied loads on a membrane geometry entirely included in the volume of the membrane. To match this objective, the algorithm starts from a parametric description of the stress potential and follows two subsequent optimization steps: in the first one, the parameters of the stress potential are optimized, through the solution of the transverse equilibrium equation of a membrane, assuming fixed boundary conditions for the membrane shape; in the second step, the boundary values of the membrane shape are selected, by minimizing an error function defined as a quadratic approximation of the characteristic function over the planform of the membrane.
The main advantage of this approach is the possibility of solving, in each optimization phase, an elliptic problem, since the stress potential is assumed to be known and concave. This strategy permits the solution of problems for which a concave shape, comprised in the volume of the vault, does not exist, such as in the case of cross vaults. For such vaults, in fact, the search of a stress field starting from a possible membrane shape (in-between the intrados and the extrados of the vault) is described by a non-elliptic second order partial differential equation, that cannot be solved by the simple inversion of the matrix arising from its discretization.
The algorithm is applied to two benchmarks: a cloister vault, with the dimensions of a vault in Palazzo Caracciolo di Avellino, in Naples, and the cross vault of the anterior portico of San Pietro in Vineis, in Anagni (Italy). For both problems, given a suitable parametric description for the stress potential, a membrane shape is found, fulfilling the requirements of Limit Analysis with acceptable geometrical safety factors, demonstrating the validity of the proposed methodology in assessing the stability of vaults.
Future perspectives for this work will include the development of an algorithm for which the stress potential and the membrane shape are optimized simultaneously, under the constraints of the equilibrium and of the Safe Theorem of Limit Analysis.