Optimal archgrids: a variational setting

The paper deals with the variational setting of the optimal archgrid construction. The archgrids, discovered by William Prager and George Rozvany in 1970s, are viewed here as tension-free and bending-free, uniformly stressed grid-shells forming vaults unevenly supported along the closed contour of the basis domain. The optimal archgrids are characterized by the least volume. The optimization problem of volume minimization is reduced to the pair of two auxiliary mutually dual problems, having mathematical structure similar to that known from the theory of optimal layout: the integrand of the auxiliary minimization problem is of linear growth, while the auxiliary maximization problem involves test functions subjected to mean-square slope conditions. The noted features of the variational setting governs the main properties of the archgrid shapes: they are vaults over a subregion of the basis domain being the effective domain of the minimizer of the auxiliary problem. Thus, the method is capable of cutting out the material domain from the design domain; this process is built in within the theory. Moreover, the present paper puts forward new methods of numerical construction of optimal archgrids and discusses their applicability ranges.


Introduction
Funiculars are planar frameworks which, albeit subject to a transverse and transmissible load, do not undergo flexure, i.e. they remain bending-free and, consequently, are not subject to transverse shear forces, while the axial stress resultants are characterized by a fixed sign: all the bars are in tension (or-in compression). The load is called transmissible, if it follows the design, yet keeping the given (usually vertical) direction and keeping the value of its intensity q(x) measured per unit length of the line orthogonal to the direction of the load [q] = [N/m], cf. Fuchs and Moses (2000). Of the funiculars being fully and uniformly stressed (up to a given limit (equal, say, −σ C , σ C being the permissible stress in compression)) one can find one of the least volume. Its rise cannot be too high but also cannot be too small, since big values of thrusts (i.e. horizontal reactions at the pin supports) lead to an inevitable increase of areas of the cross-sections, resulting in the increase of the volume. Thus, the rise should be appropriately chosen. Let z = z(x) represent the shape of the funicular pin supported at A and B at the levels z(0) = h A , z(l) = h B .
The optimal rise of the least-volume funicular is determined by the condition 1 l l 0 dz dx discovered by Rozvany and Prager (1979) and Rozvany and Wang (1983) and called there: the mean square slope condition. An illustrative example of the optimal funicular corresponding to a particular transverse load is shown in Fig. 1. The level function of a funicular follows the diagram of the bending moment in a simply supported beam with supports of the same ordinates as points A and B and subject to the same load q(x). Thus, for a given vertical transmissible load of intensity q(x), referred to the unit area of the basis plane, one can construct a family of funiculars and select one of the least volume, see Section 6.1 in Lewiński et al. (2019a). The main conclusion is that in a given plane one can construct frameworks, usually having Fig. 1 The funicular of optimal shape corresponding to the given load. Its level function satisfies the condition (1) forms of arches, which are simultaneously tension-free and do not undergo bending. Admitting tension leads to Michelllike designs of great complexity, see Darwich et al. (2010). A similar theoretical construction in space is much less obvious. A natural counterpart of a planar bending-free arch is a membrane shell. Of three local equilibrium equations, one is algebraic: where N 1 , N 2 represent the membrane stress resultants corresponding to the principal curvature parameterization, while the principal curvatures are denoted by R 1 , R 2 ; q represents intensity of normal load, [q] = [N/m 2 ]. Thus, if the normal load is absent and the curvatures are positive, the stress resultants N 1 , N 2 must have opposite signs. At least in this case, a membrane shell cannot be subject only to compression and work in a bending-free state. On the other hand, we remember the solutions of statics of membrane shells of revolution under self-weight: at the pole, both the membrane stress resultants are negative and equal but along the meridian, the diagrams of these stress resultants have different distribution: the meridional force (i.e. its absolute value) increases, while the initially negative circumferential stress resultant increases and, crossing a zero value at a certain circumference, attains positive values along the support, see Flügge (1960). These examples show that, in general, one cannot construct a membrane shell which under a transmissible vertical load would work in a tension-free state. That is why instead of membrane shells, the tensionfree grid-shells are preferred in many applications. The grid-shells are tension-free spatial frameworks formed on a single or on several surfaces, their form being selected such that the bending is vastly eliminated, see Richardson et al. (2013) and Jiang et al. (2018). Grid-shells are usually constructed in two steps. First, the geometry of the deformed configuration is found, e.g. by dynamic relaxation, see Day (1965). Then the optimal position of nodes is found to minimize the total weight and impose additional design conditions, see Richardson et al. (2013). The recent paper by Jiang et al. (2018) discusses three more efficient methods of constructing grid-shells: the force density method (FDM) by Linkwitz and Schek (1971) and Schek (1974); the potential energy method (PEM); and the ground structure method (GSM). In the FDM, the values of the member forces are predefined, while the nodal position is determined by the equilibrium conditions in the undeformed state. Thus, the designer has a full control over the values of axial stresses and it is possible to eliminate tension (or compression) a priori. The GSM is also capable of eliminating tension (or compression) but the solution does not constitute a single surface structure; in general, a multi-storey structure appears, cf. Lewiński et al. (2019b).
An alternative method of construction of grid-shells can be inferred from the theory of archgrids, proposed by Rozvany and Prager (1979) and extended in Rozvany et al. (1980Rozvany et al. ( , 1982, Rozvany and Wang (1983), and Wang and Rozvany (1983). This theory has been recently revisited in Lewiński et al. (2019a, Sec. 6.3.6;2019b), Czubacki and Lewiński (2019), and Czubacki et al. (2018). Archgrids are composed of infinitely thin arches formed in the mutually orthogonal planes: x = const., y = const., orthogonal to the basis Ω corresponding to z = 0. The given load of intensity q is transmissible along the z axis and is decomposed into the loads q x , q y acting on the arches in y = const. and x = const. planes. The shape of the arches is chosen such that the bending is fully eliminated. Consequently, also the transverse shear forces vanish. Thus, the state of stress in the arches is axial: the only stress components are normal to the cross-sections. One can design the areas of cross-sections of the arches such that the state of stress is fully uniform within the arches and equals −σ C . Of such uniformly stressed structures, one can find the one of the smallest volume. Rozvany and Prager (1979) discovered that this least volume structure is formed on a single surface. More precisely, if z = z x (x, y) represents the surface composed of arches formed in y = const. planes, and if z = z y (x, y) represents the surface composed of arches formed in x = const. planes, then the conditions of minimum volume imply the identity: z x (x, y) = z y (x, y) = z(x, y) there, where both the arches are subject to non-zero loads. On the other hand, arches in one direction can also appear and even some part of the basis Ω can be uncovered. The point loads are also allowed and lead to arches of finite cross-sections. Nevertheless, the optimal structure is a grid shell in all cases; the problem of its construction belongs to the class of topology optimization problems.
The numerical solutions shown in the mentioned papers by Rozvany and co-authors have been constructed iteratively, the optimization problem itself had not been explicitly formulated. Czubacki and Lewiński (2019) noted that the archgrid problem reduces to the two mutually dual problems, to some extent similar to those constituting the scalar Beckmann problem, see Santambrogio (2015), forming the Michell's problem, see Lewiński et al. (2019a) as well as the free material design problems, see Bouchitté and Buttazzo (2001), Czubacki and Lewiński (2015), and Czarnecki and Lewiński (2017). The minimization problem is characterized by the functional of linear growth; its dual version is a maximization problem of the virtual work with locking conditions. Just this decomposition governs the topology optimization features of the setting and of the final form of solutions. Moreover, the numerical and mathematical methods to attack the problem should be appropriately chosen, as for functionals with integrands of linear growth.
The mathematical theory of archgrids evenly supported along the closed contour of the basis domain has been put forward in Czubacki and Lewiński (2019). The present paper is aimed at extensions of this theory towards archgrids unevenly supported. The level of the edge support is determined by the values of the ruled surface given by the equation: z 0 (x, y) = a 1 xy + a 2 x + a 3 y + a 4 . In the first step, the surfaces z x (x, y), z y (x, y) are generated by the family of planar arches in the sections y = const., x = const., supported along the contour Γ of the domain Ω on the level z 0 along Γ . These surfaces determine the grid shells carrying the loads q x (x, y), q y (x, y) such that q(x, y) = q x (x, y) + q y (x, y). The decomposition of the given load q into q x , q y is one of the unknowns of the archgrid theory. Upon introducing the areas of crosssections for assuring a uniform state of stress in the whole structure, we state the problem of the volume minimization. This leads to various optimality conditions, one of them being the identity z x (x, y) = z y (x, y) = z(x, y) valid in these subdomains, where two families of arches carry the load. All the optimality criteria will be discussed and the mathematical formulation of the optimum design problem will be stated; it turns out that the level function z(x, y) of the vault is governed by a well-posed problem composing of two auxiliary mutually dual problems. Upon finding the static unknowns, one can compute the optimal areas of the cross-sections to assure the uniform distribution of stress. The optimum archgrid is the lightest among all statically admissible and fully stressed.

Construction of the bending-free arches
The present section shows that for any transverse and transmissible load, one can construct a three-hinge arch which is bending-free.
Consider a three-hinge planar arch subject to a vertical load of intensity q(x), 0 ≤ x ≤ l, referred to its projection on the horizontal axis x. The hinges at the support lie on the levels: z = h A and z = h B , see  Along with the arch problem, an auxiliary beam supported at points lying below the supports A and B of the arch, of ordinates x = 0 and x = l, is considered. The beam is subject to the same load q(x), 0 ≤ x ≤ l.
The vertical reaction at the left support A equals The bending moment in this beam at the section of ordinate x is expressed by The condition of the total equilibrium of the arch reads: where H represents the horizontal reaction at the support A.
We consider only such loads q for which H > 0. Note that the vertical reactions of the arch and the auxiliary beam are linked by where tan(α) = (h B − h A )/ l. The bending moment in the arch at the section of the ordinate x can be expressed by The shape of the arch can now be adjusted to the given load by requiring that M(x) = 0 for 0 ≤ x ≤ l. This condition determines the shape of the bending-free arch: Thus, where T represents the transverse shear force in the auxiliary beam of Fig. 2b. This transverse shear force satisfies the differential equilibrium equation Its variational counterpart has the form where V 0 represents the set of trial functions vanishing at x = 0 and x = l; the regularity conditions involved in the definition of V 0 will not be discussed in the present paper. The equalities (9) and (10) imply called the funicular equation.
One can prove that the arch of the form given by the elevation function (8) is not only bending-free but is also not subject to transverse shear. It is subject only to compression. The axial force N(x) in this arch is linked with the horizontal reaction H by the formula where ϕ(x) is the angle of inclination of the tangent at point (x, z(x)) to the axis x, hence tan(ϕ(x)) = dz/dx, cf. Fig. 3. Thus, N(x) < 0 along the whole arch. The areas A(x) of the cross-sections can be chosen such that the normal stresses are made uniform within the whole arch and equal −σ C , σ C being the permissible stress in compression. To this end, we assume |N(x)|/A(x) = σ C , or

The bending-free arch of least volume
Among the arches which are bending-free and uniformly stressed, one can find an arch of the least volume. The volume of the arch of areas of cross-sections given by (14) equals By using the formulas 1 cos 2 (ϕ(x)) = 1 + tan 2 (ϕ(x)), and considering that the integral of T (x) along the span of the auxiliary beam vanishes, we obtain wherel = (1 + tan 2 (α))l and · stands for the norm 1 The optimum design problem (P 0 ) Among bending-free three-hinge arches corresponding to a given transmissible load q(x) find the least-volume arch.
The solution.
All such arches are characterized by the level functions z(x) given by (8), involving the design variable H . By minimizing the volume V over H , we obtain the optimal value of the horizontal reaction H =Ȟ , oř while the volume of the lightest arch equalš Its elevation function z =ž is given by (8) and (19), oř see Fig. 1. Let us differentiate both sides of (21): The right hand side has the unit norm; consequently, the norm of the left hand side is equal to 1 and this condition leads to, cf. (1) This condition means that the optimal arch cannot be too shallow and also cannot be too high. If it were very shallow, the load would cause great values of the thrust H ; consequently, the areas of the cross-sections would be appropriately big to assure a constant normal stress. On the other hand, if an arch is higher and higher than its volume increases. The condition (23) indicates that the rise of the arch should be perfectly chosen. It is worth noting that this is an integral condition and not a point-wise one. 2 Moreover, note that the bound increases along with the increase of α. Thus, if the difference between the levels of the supports is bigger and bigger, the mean rise increases accordingly. This fact is not obvious and can serve as an hint for a designer. The arch of minimal volume has the shape given by the elevation function (21), while the areas of its crosssections are given by (14) and (19). The optimal areas of cross-sections are expressed as below Note thatǍ cannot degenerate along the arch and assumes finite values.
Remark 1 It is worth noting that the optimal volumeV (given by (20)) is proportional to the potential energy of the load, namely The proof runs as follows. Sincez vanishes at the ends of the beam, we can choose v =z in (11), which implies Substitution of (27) into (26) gives which ends the proof.
Remark 2 Consider the following two auxiliary problems of the variational calculus. Problem (P 0 ) Find the minimizer Q =Q of the problem Problem (P 0 ) Find the maximizer v =v ∈ V 0 of the problem If regularity assumptions are appropriately chosen, one can prove that Z 1 = Z 2 , cf. Czubacki and Lewiński (2019). Moreover,v andQ are linked by andQ represents just the transverse shear force appearing in the auxiliary beam of Fig. 2b, orQ = T . Thus, where M is bending moment in the auxiliary beam of Fig. 2b. Moreover, according to (21), the elevation function of the optimal arch can be expressed by which impliesz(x) = lv (x). It is seen that the fielď v(x) attains the bound in the condition involved in (30): dv/dx = 1.
In the problem considered, the variational problems (P 0 ), (P 0 ) are not helpful in construction of the solution, since the solution (21) has been directly constructed. Nevertheless, the formulations of problems (P 0 ), (P 0 ) are inspiration for the construction of the optimal archgrids which are the main subject of the present study.

Archgrids unevenly supported along a closed contour
The aim of this section is to put forward a variationally consistent theory of optimal archgrids.

The stress-based construction of archgrids
The term archgrid is indissolubly bonded with its leastvolume property. It will be shown that just the least-volume condition leads to shaping a single surface on which the arches (in the planes x = const., y = const.) are formed. The basis of the archgrids to be designed will be a plane domain Ω parameterized by a Cartesian coordinates system (x, y). Assume that the straight lines y = const. cut the contour Γ of the domain Ω at two points: Fig. 4. Assume that the domain Ω is included in a rectangle of dimensions a 0 and b 0 and let the point (0, 0) be one of its vertices; The sides of the rectangle are given by Fig. 4. Let us form over a given rectangle a ruled surface whose elevation function is expressed by The angles of inclination of the lines x = const., y = const. are expressed by: The vault will be designed over the basis Ω, along the contour Γ . Its elevation along the contour Γ will be given: The elevations of the vault at the points of the contour Γ : The following equalities hold for the subsequent cross-sections y = const., x = const., see Fig. 5. The vault to be designed is assumed as an archgrid composed of two families of arches in the directions y = const., x = const., hence formed in the mutually orthogonal planes. The given load q, referred to a unit area of a basis plane, The arches in the sections y = const., subjected to the vertical load of intensity q x (x, y), are formed such that they are bending-free with using the method of Section 2.

Fig. 5
The ruled surface z = z 0 (x, y) over the rectangle surrounding the basis domain Ω Consider a simply supported beam of length l x (y). The load q x (x, y) causes the transverse shear force T x (x, y) and the bending moment M x (x, y). These fields satisfy the differential equations of equilibrium: The level function of the arch is assumed such that the arch is bending-free. Thus, according to (8), the ordinates of the arch formed for y = const. are given by In the cross-sections x = const., the optimal arches, subjected to the vertical load of intensity q y (x, y), are constructed in a similar manner. Consider now simply supported beams of lengths l y (x). The load of intensity q y (x, y) causes the fields T y (x, y), M y (x, y) satisfying The elevation function of the arch is assumed by the rule (8), or The arch thus constructed is bending-free. The horizontal reactions H x (y), H y (x) are assumed to be non-negative. Note that the fields T x , T y satisfy the equilibrium equation At this stage of the study, the designed vault is formed on two surfaces of elevation functions: z x (x, y), z y (x, y). Their values coincide along the contour of Γ .
The values of z 0 along the contour are determined by the ruled surface (34). In the sequel, our aim is to construct a vault of the minimal volume. It will occur that this condition will result in the identity z x = z y within Ω. Just this identity will be interpreted as a junction of both the families of arches to form a one shell-like structure.
The computation of volumes of the both families of arches requires analysis of their slopes. By differentiating (40) and (42), we get Let ϕ x , ϕ y represent the inclination of the tangents to the arches in the planes y = const., x = const., respectively, see Fig. 6. Thus, To make the distribution of the axial stress uniform and equal −σ C , we shall choose the densities of areas of crosssections in an appropriate manner. By using the rule (14), one obtains The elementary volumes of the arches are expressed by where ds x = dx cos(ϕ x ) , ds y = dy cos(ϕ y ) Substitution of (47), (49) into (48) gives The total volume of the arches parallel to the axis x is expressed by while the total volume of the arches parallel to the axis y equals Substitution of (46), (45) into (51) results in which reduces to wherẽ l x (y) = l x (y) 1 + tan 2 (α(y)) In a similar manner, we obtain wherẽ l y (x) = l y (x) 1 + tan 2 (β(x)) (57) Fig. 6 Construction of the archgrid from two families of arches in the planes y = const., x = const. . The auxiliary beams are subject to the loads q x , q y , respectively Now we are ready to formulate the main problem of the present paper, namely the problem of constructing the leastvolume archgrid: The minimizing pairs (Ȟ x ,Ȟ y ), (Q x ,Q y ) are stationary points of the Lagrangian L = L(H x , H y , T x , T y , λ) given by (58) where λ(x, y) is a Lagrangian multiplier corresponding to the equilibrium condition (43) and vanishing along the contour Γ . By integrating by parts, we rearrange the Lagrangian to the form where Let us note that the condition is equivalent to (43). Thus, the fields T x , T y involved in (59) will be treated as independent. Moreover, v = λ/2 plays in (59) the role of a Lagrangian multiplier. We require δL = 0 with respect to independent variations δv, δH x , δH y , δ T x , δ T y . In this way, the stationary conditions of the functional (59) are obtained: the condition (61) and where the norms ρ x (· ; ·, ·), ρ y (· ; ·, ·) are defined by The stationarity conditions (62) and (63)  If H x , H y are given and H x > 0, H y > 0, then the solution v of problem (P 1 ) is unique. This function v determines uniquely the fields T x , T y by (65). The formulae (62) may be re-written as where now v is the solution to the problem (P 1 ). By combining (68) with (45), one gets (69) hence, The functions s 1 (y), s 2 (x) can be chosen such that z x = z y in the whole basis domain including the boundary. Indeed by taking one obtains the identity (tan(α(y)))x + s 1 (y) = (tan(β(x)))y + s 2 (x) cf. (36). Since v vanishes on Γ , we note that z = z x = z y assumes the values z 0 (x, y) for (x, y) lying on Γ . Thus, if H x (y), H y (x) are given, then the function determines the shape of the vault. To assure that this shape refers to the structure of the minimal volume, the functions H x (y), H y (x) should be appropriately selected. The relations (63) make it possible to express the total volume V = V x + V y in terms of the fields T x , T y . The direct substitution of (63) into (54), (56) results in where for w = (u, v) we define The vector field T = ( T x , T y ) is subject only to the condition (43) (or equivalently: (61)). Now the problem (P) of minimization of the volume of the archgrid can be explicitly formulated as (cf. the theoretical scheme given in Fig. 7): The problem above is well posed since the set of vector fields Q of given divergence is non-empty. Moreover, let us note that the functional ℘ (w) possesses the properties of a norm of a vector function of argument w = (u(x, y), v(x, y)).
Since the functional in (76) is convex and the set of admissible vector fields Q is non-empty, the problem (P ) is well-posed. The problem of existence depends on regularity of the data.
LetQ Q Q = (Q x ,Q y ) be the minimizer of (76). By substituting T =Q Q Q into (63), one obtains the optimal distributionsȞ x (y),Ȟ y (x) of the boundary forces. Substitution of these results into (67) fixes the coefficients in the bilinear form aȞ (·, ·) governing the optimal solution. IfȞ x > 0,Ȟ y > 0, then this bilinear form is positive definite and the problem (P 1 ) possesses a unique solutionv which determines the surface of the archgrid by (73).

The displacement-based construction of the archgrid
It turns out that the stress-based problem (P ) has its displacement-based counterpart and just the latter problem will be, usually, easier to solve. The aim of this section is to put forward this dual formulation.
Let us start by noting that the gradient ofv can be computed directly from (62) by putting Let us substitute: H x =Ȟ x , H y =Ȟ y , and into (63). IfȞ x > 0,Ȟ y > 0, then the relations (63) imply The fieldv determined by the solution of the problem (P 1 ) with the optimal values of the coefficients H x =Ȟ x , H y = H y of the quadratic form (67) is the maximizer of Indeed, the problems (P ), (P ) are mutually dual; Z 1 = Z 2 , or the duality gap vanishes. For the proof, the reader is referred to Czubacki and Lewiński (2019). The proof applies here since the formal structure of problem (P ) and problem (78) in Czubacki and Lewiński (2019) is the same. It is sufficient to replace l x , l y withl x ,l y and if Q Q Q = (Q x ,Q y ) is the minimizer of (76) and ifv is the maximizer of (80), then these fields are linked by the optimality conditions HavingQ x ,Q y , one can find ∇v but not vice versa; the relations (81) are non-invertible. Therefore, the construction of the optimal archgrid necessitates solving the two independent problems: (P ), (P ), while (81) can play the role of a check. Upon findingv, one can construct the level functionž of the optimal archgrid by the rulě where z 0 (x, y) is given by (34). One can prove that the functionž(x, y) satisfies the equalities: x 2 (y) These equations are counterparts of the condition (23) known from the theory of planar funiculars. Let us prove the first of the conditions (83). To this end, we insert (cf. (69)) into the first of relations (79). One finds While computing the l.h.s of (85), one can make use of the formulas: z(x 2 (y), y) −ž(x 1 (y), y) = l x (y) tan(α(y)) and thus rearrange (85) to the form of the first of equalities (83). One can say that the conditions (83) introduce the upper bounds on the slope of the arches in both directions. These bounds are increasing functions of arguments α(y), β(x). The theoretical construction of the optimal archgrids governed by problem (P) is reduced to solving problems (P ), (P ), their juxtaposition being presented in the diagram in Fig. 7.
Having solved the problems (P ), (P ), one can determine the areas of arches' cross-sections, which will be discussed in the sequel.

Determination of the densities of the areas of arches' cross-sections
The solution to the problem (80) determines the shape of the optimal archgrid according to (82) but, in general, does not allow for a direct computation of the densities of areas of cross-sectionsǍ x ,Ǎ y in the optimal structure. These fields are determined by the solution (Q x ,Q y ) to the problem (76). To find the explicit formulas for the densities of the Fig. 7 The tabular scheme of the theory of optimal archgrids area cross-sectionsǍ x ,Ǎ y , we shall make use of (47). The forces H x , H y involved in these formulas are found by putting T x =Q x , T y =Q y in (63).
Having found the minimizerQ Q Q = (Q x ,Q y ) of problem (P ) and the maximizerv of problem (P ), one can compute the volume of the lightest archgrid either by or by which shows that the optimal volume is proportional to the potential energy of the loading.

Optimal archgrids over rectangular domains-the complete numerical approach with using Legendre polynomials
The present section is aimed at numerical constructions of optimal gridworks over rectangular bases, corresponding to selected transmissible loads including the loads concentrated along curves. The construction is based on tackling simultaneously the problems (P ) and (P ) along with checking the optimality conditions (81). Both the problems (P ), (P ) are solved numerically by using Legendre polynomials.

Construction of the solution to the problem (P )
We confine consideration to the case of the bases being rectangular domains The nth Legendre polynomial is given by where [·] is a function which rounds down its argument to the nearest integer. Let us introduce polynomials T n (t), called anti-derivatives of P n (t), cf. Abedian and Düster (2019), such that: and T n (−1) = T n (1) = 0 for n ≥ 1. The explicit form of T n reads while T 0 (t) = t does not satisfy the above boundary conditions. The test functions in (80) will be expressed by the truncated series involving the T n functions, as below withx = x/L x ,ỹ = y/L y ; v ij are unknown coefficients, the numbers m x , m y control accuracy of the method. Due to (94) the partial derivatives of v assume the form: hence, . 8 The numerical flowchart. The unshaded (blank) blocks refer to the steps performed by symbolic computations, while the shaded blocks refer to the steps programmed numerically or ρ y ∂v ∂x ; x 1 (y), x 2 (y) and the expression ρ x ∂v ∂y ; y 1 (x), y 2 (x) can be written in a similar manner. Note that a significant number of components at the r.h.s. of (99) cancel by orthogonality of Legendre polynomials: where δ ij denotes the Kronecker delta. This simplifies the formulation (80) essentially and paves the way for an efficient programming of the maximization problem, see the r.h.s. of Fig. 8. The boundary condition v = 0 on Γ is satisfied identically. The coefficients v ij are the only unknowns of the problem. The constraints in (80) are satisfied along the lines: x = x i , i = 1, . . . , s x and along the lines y = y j , j = 1, . . . , s y , where the numbers s x and s y are fixed and chosen simultaneously with numbers m x , m y . For each x i , i = 1, . . . , s x , the corresponding slopes β(x i ) and the auxiliary lengthsl y (x i ) are computed by (37) 2 and (57) respectively. For each y j , j = 1, . . . , s y , the corresponding slopes α(y j ) and the auxiliary lengthsl x (y j ) are computed by (37) 1 and (55) respectively. In the problems considered in the present paper, the choice: m x = m y = m occurred to be efficient.
The numerical approximations to the definite integrals are performed with using the NIntegrate function available in Wolfram Mathematica software.

Construction of the solution of the problem (P )
The unknowns of the problem (P ) are expressed by using the Legendre polynomials P n and their anti-derivatives T n as follows and the vertical load is represented by where (103) The differential condition nested in (76) makes it possible to eliminate the unknowns Q y ij by the algebraic equation thus making the unknowns Q x ij the only design variables of the problem. The subsequent steps of the method are schematically shown in the flowchart in Fig. 8. Upon finding the minimizer Q x ij , Q y ij , one can either compute the volume by (91) or multiplying by the factor 2 the objective function F , see l.h.s. of Fig. 8. This method delivers directly the distribution of the horizontal forces along the edges, see (63). However the method does not provide a direct formula for the optimal level functionž(x, y). To find it, one should integrate the equalities (62) and augment them with the boundary conditions. Let us emphasize that the optimal shape is directly delivered by solving problem (P ).
For the same case of λ = 0, the results obtained by the method (P ), for selected values of n x = n y = n, are listed in Table 2. The values of volumes obtained by the method (P ) corresponding to the increasing number m x = m y = m form an increasing sequence (a lower converging sequence), see the first row in Table 1, while the sequence of volumes obtained by the method (P ) for increasing number n x = n y = n is a decreasing sequence (an upper converging sequence), see Table 2. Both sequences tend to the same optimal value, thus confirming the duality gap between problems (P ) and (P ) being zero.
Of two methods shown in Fig. 8, the method based on (P ) produces the sequence of results converging faster. That is why the results in Table 2 could be confined to n = 30, which corresponds to the accuracy of the method (P ) for m = 75.
The method (P ) gives the optimal shapes of optimal archgrids directly. They are not concave, cf. Fig. 9. Let us emphasize thatž (λ =0) = z 0(λ) +v (λ=0) . On the other hand, the method (P ) delivers directly the distribution of the horizontal reactions H x , H y and the areas densities A x , A y . The method (P ) for m = 75 gives the results of the gradient of v of better accuracy than the accuracy with which the method (P ) for n = 30 delivers the quotients Q x /H x , Q y /H y , see Fig. 10 corresponding to the case λ = 0. The same figure shows the precision with which the optimality conditions (81) are fulfilled. To measure the relevant error, the function:  is introduced and the function error y (x) is defined similarly. These errors, corresponding to the sections x = const., y = const., depend on the distance to the edges. The errors are small far from the edges and increase in the boundary zone, see Table 3. This phenomenon appears because in the problem (P ), the boundary conditions are not explicitly imposed.

The optimal archgrids-other illustrative designs
The optimal archgrids to be constructed in this section will be found by solving the problem (P ) in order to have directly the optimal level functionž(x, y). Two kinds of vertical loads will be discussed: the load concentrated along

Fig. 11
The square domain 2L × 2L as a basis and the contour of a circle of radius R = 0.75L along which the load p is applied (a). The archgrid is supported such that z = 0 along Γ (case i), or z = z 0 along Γ , z 0 being shown in (b) (case ii)  a contour of a circle and the load uniformly distributed within a circle. Fig. 13 The cross-sections of the level functionž of the optimal archgrid shown in Fig. 12 for y = i 25 L, i = 0, . . . , 24. The red dashed lines correspond to non-optimal arches for which the bounds involved in (80) are not saturated 5.1 The optimal archgrids over a square domain for the load uniformly distributed along a contour of a circle The optimal archgrids will be constructed over a square domain Ω, of side of length 2L, loaded uniformly along a contour of a circle of origin in the middle of domain and a radius R, see Fig. 11. The load acts along the z axis orthogonal to Ω, the intensity of the load being p = const., [p] = [N/m]. The two cases of the support will be considered: (34) with (c 1 , c 2 , c 3 , c 4 ) = (L, 0, L, 0), see Fig. 11.
Due to symmetries of the problem v kl = 0 for k, l being even numbers and v kl = v lk due to symmetry of a square domain. Thus, integral constraints in (80) are applied only for x i = (i − 1) L s for i = 1, . . . , s and s = 2m. The accuracy of the result of the optimal volume is controlled by the number m. For increasing m, the volume tends to cf. Table 4. The shape of the optimal archgrid supported as in case (i) is surprising, see Fig. 12a. The vault is formed not over the whole basis domain; the four squares of the basis: Fig. 12d. The load acting along the circle is transmitted to the sides of the basis domain by straight arches forming four ruled surfaces over the domains: Fig. 13. The net of bars over the Fig. 14 The contour plot of the optimal solution of problem (i), Section 5.1. The given load of intensity p is applied along the circular contour shown as a red dot-dashed line. The blue dashed lines indicate the directions of the optimal arches Fig. 15 The shape of the optimal archgrid for the problem (ii) of Section 5.1. The axonometric view interior square B 1 B 2 B 3 B 4 forms a surface of non-zero Gauss curvature, see Fig. 14. This surface is continuously connected with the ruled surfaces mentioned above by the ruled surfaces over the domains: B 1 C 1 B 2 , . . . , B 4 C 4 B 1 . Note that although the archgrid is not concave, all the arches are in compression and the stress distribution is uniform, as assured by appropriate design of areas densities (89). Let us stress that the optimization process has determined the position of the load p over the contour C 1 C 2 C 3 C 4 . The points B 1 , . . . , B 4 are at the highest levels (see Fig. 14), which may suggest that the load is concentrated at these vertices, which is obviously not the case. Let us stress here that the exact analytical solution of this result is still pending. Fig. 16 The cross-sections of the level functionž of the optimal archgrid shown in Fig. 15 for y = i 25 L, i = 0, . . . , 24. The red dashed lines correspond to non-optimal arches for which the bounds involved in (80) are not saturated Fig. 17 The contour plot of the optimal solution of problem (ii), Section 5.1. The given load of intensity p is applied along the circular contour shown as a red dot-dashed line. The blue dashed lines indicate the directions of the optimal arches The optimal archgrid for the case (ii) of the support has two symmetry planes, see Figs. 15 and 16, where the crosssections y = const. are shown. We note that the arches touching the contour of the basis form the ruled surfaces as in case (i). These ruled surfaces and the surface over the interior region B 1 B 2 B 3 B 4 are connected by other four ruled surfaces, see Fig. 17. The optimal volume equals V /V p = 5.361, V p = pL 2 /σ C (107) and this value is bigger than that of case (i) by 1.26 %. Thus, the construction is heavier, as expected.

The optimal archgrid over a square domain for the load uniformly distributed within a circle
The aim is to construct an optimal archgrid over a square domain of the side of length 2L for the load of intensity q uniformly distributed within a circle of origin at (0, 0) and radius R = 0.75L, see Fig. 18a. The level functionž is constructed by using Legendre polynomials by the method described above. The optimal shape of the archgrid (Fig. 18b) is similar to that of Fig. 12a, but now the central part is concave and not convex as before. This structure is composed of four ruled surfaces A 1 A 2 B 2 C 1 B 1 , . . . , A 7 A 8 B 1 C 4 B 4 extended towards the interior of the basis domain by the surfaces B 1 C 1 B 2 , . . . , B 4 C 4 B 1 -which here are not ruled surfaces-and completed by the surface of non-zero Gauss curvature over the B 1 B 2 B 3 B 4 square (Fig. 19).

Final remarks
The theory of optimum design of trusses teaches us that the condition of minimum compliance with a prescribed value of the volume leads to designs in which all the bars are stressed up to the same given limit (cf. Lewiński et al. (2019a), Section 2.1). This suggests that the design optimization of archgrids based on the uniform stress condition should lead to the designs of least compliance. Yet this conjecture needs deeper studies due to issues concerning deformation behaviour of archgrids; note that a single funicular is geometrically variable, which makes the kinematic analysis nontrivial. The identity of both the surfaces z x (x, y) = z y (x, y) of two families of arches shown in Section 3 refers obviously to the undeformed state. Yet the arches of both the families are not connected at nodes on this surface. Consequently, they may form two different surfaces upon deformation. Yet the choice of the areas of the cross-sections of both the Fig. 18 The square domain 2L × 2L as a basis and the circular domain of radius R = 0.75L where the uniform load of intensity q is applied (a). The axonometric view of the optimal archgrid (b). The crosssections of the level functionsž for y = i 25 L, i = 0, . . . , 24. The red dashed lines correspond to non-optimal arches for which the bounds involved in (80) are not saturated (c). The upper view of the layout of arches (d) Fig. 19 The contour plot of the optimal solution of the problem shown in Fig. 18. The given load of intensity q is applied within the circular contour shown as a red dot-dashed line. The blue dashed lines indicate the directions of the optimal arches families of arches assures uniform stress distribution. We conjecture that this property singles out the rational design which assures the required deformation compatibility. In the forthcoming papers, this question will be addressed by a rationalization of archgrids by replacing them by frames of finite number of arches connected by stiff joints. Like Michell structures (see He and Gilbert (2015)), the archgrids may be effectively rationalized, thus showing their lightness, stiffness, and applicability in the civil engineering practice.
The theory put forward concerns the archgrids composed of two families of arches in x = const and y = const planes. Such constructions can be improved, noting that the load applied close to the corners is better transmitted by arches in the planes going in ±45 • directions, cf. Sec.6.3 in Lewiński et al. (2019a) and examples in Lewiński et al. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creativecommonshorg/licenses/by/4.0/.