Sparse Stress Structures from Optimal Geometric Measures

Identifying optimal structural designs given loads and constraints is a primary challenge in topology optimization and shape optimization. We propose a novel approach to this problem by finding a minimal tensegrity structure—a network of cables and struts in equilibrium with a given loading force. Through the application of geometric measure theory and compressive sensing techniques, we show that this seemingly difficult graph-theoretic problem can be reduced to a numerically tractable continuous optimization problem. With a light-weight iterative algorithm involving only Fast Fourier Transforms and local algebraic computations, we can generate sparse supporting structures featuring detailed branches, arches, and reinforcement structures that respect the prescribed loading forces and obstacles.

Fig. 1.Given loading forces f  at some nodes p  in the presence of obstacles, find a minimal graph of struts and cables that form a structure capable of supporting the forces while avoiding obstacles.
Designing the geometry for a sparse support structure is an important task in topology optimization and structural design, yet it remains a computationally challenging problem.Typically, the problem is formulated as a large-scale parameter optimization problem over an elastostatic solver or as a combinatorial problem [Rozvany 2009].In this paper, we approach the subject matter by considering a geometric graph optimization problem: Find a minimal network of curves connecting a given set of ends while avoiding a given set of obstacles, so that it represents a minimal supporting structure for a given set of loading forces given at the ends.
We refer to this problem as the minimal stress reconstruction problem.Using geometric measure theory, we demonstrate that this problem can be formulated as a simple continuous optimization problem.We also show that approximate solutions to this problem exhibit intricate emergent geometric structures and patterns such as branches, trusses, and arches (Fig. 2).Despite the simplicity of the mathematical model, the phenomenological richness of its results suggests numerous potential applications in topology optimization, architecture and tensegrity structure design, as well as theoretical connections to the study of branched optimal transport [Xia 2010;Bonafini et al. 2018;Pegon et al. 2019] and cytoskeletal networks [Stamenović et al. 1996;Ingber 2003].This paper focuses on the minimal stress reconstruction problem in the plane with arbitrary user-prescribed obstacles and weights.The optimization algorithm is a simple iterative scheme involving only fast Fourier transforms and local calculations.

RELATED WORK
Topology Optimization.Topology optimization (TO) aims to find topology and shape for a structure using a given volume of material which also minimizes a physical energy (typically stress-strain energy).The diverse range of methods used for TO constitute a wide and deep body of literature [Bendsoe and Sigmund 2003;Rozvany 2009;Sigmund and Maute 2013;Deaton and Grandhi 2014], and include greedy heuristic methods, levelset optimization methods, and genetic algorithms [Rajeev and Krishnamoorthy 1997].TO typically manages the material volume using explicit constrains on the total density.Our method instead directly minimizes the support of our stress tensor, and thus describes a different class of optimization problems.Compared to traditional TO, our model is simpler and depends on fewer parameters.
Our method is closely related to a variant of the topology optimization called the minimal stress problem, in which one minimizes the 1 norm of the stress field subject to a divergence constraint [Strang and Kohn 1983;Bendsoe and Sigmund 2003, § 3.3.3].However, unlike the classical topology optimization which produces sparse structures, a  1 -minimal stress distribution is generally a dense continuum.In this continuum the pair of principal directions of the stress matrix trail out the grid lines of the Michell structure.These minimal Michell continua are also computed on curved surfaces for architectural [Kilian et al. 2017] and mechanical [Gil-Ureta et al. 2019] applications.Our method generalizes the  1 minimization to a non-convex continuous optimization that is able to produce sparse graph structure instead of a dense continuum.
Tensegrity.Tensegrity, which is short for "tensile integrity" [Fuller 1962], is a mechanical paradigm characterized by cables and struts supporting a structure in equilibrium.It has previously been posed as a mixed integer program [Ehara and Kanno 2010], semidefinite programs [So andYe 2006, 2007;Wang and Xu 2019] and linear programs [Kushner et al. 2021], with rich connections to the geometric realization of graphs.Tensegrity also has applications to architecture [Gomez-Jauregui 2004], cell biology [Stamenović et al. 1996;Ingber 2003], robotics [Paul et al. 2006;Graells Rovira and Mirats Tur 2009], and space structures [Tibert and Pellegrino 2003].Our algorithm's output can be interpreted as a tensegrity structure whose topology has been optimized using an Eulerian method.
Varifold Methods.Varifolds are a generalization of manifolds from geometric measure theory which have appeared in some computer scientific contexts in the past.Buet et al. [2017] use "discrete varifolds" to approximate the mean curvature of input geometry.Charon and Trouvé [2013] employ varifolds to handle registration for potentially nonorientable shapes.Our method uses a varifold formulation to find potentially nonorientable and nonmanifold surfaces matching a certain force distribution.
Computational Geometric Measure Theory.A recent method uses geometric measure theoretic formulations to find area-minimizing surfaces over 3D domains [Wang and Chern 2021]; extensions to this work used deep learning to increase resolution and generalize the method to arbitrary oriented manifolds with boundary [Palmer et al. 2022].Critically, unlike these methods, our method is able to handle nonorientable submanifolds with branch points.

THEORY
Given  0 = {p 0 , . . .p  } ⊂ R2 and f 0 , . . .f  ∈ R 2 , the minimal stress reconstruction problem asks for an R 2 -embedded undirected graph ( , ),  ⊃  0 , and tensions (   )   ∈ ∈ R such that (1) Here, f  is considered to be 0 on  \ 0 .This is the condition that the graph forms a structure in equilibrium which exerts the prescribed forces on vertices in  0 .We also seek to find the graph with minimal edge length among all graphs satisfying these properties, where edge length is denoted ℓ   = |p  − p  |: minimize   ∈ ℓ   subject to (1). (2) Overview.We approach this graph optimization problem (2) by its continuous relaxation, derived in Section 3.1.In the end of this relaxation, we optimize over the continuous space of stress distributions on R 2 , instead of directly optimizing among the graph's combinatorics.The nodal forces are described as a force distribution f : R 2 → R 2 , and the edge-wise tensions are described by a symmetric matrix field as the stress distribution 1  : R 2 → R 2×2 sym .The tension on an edge   ∈  in (1) corresponds to stress being the rank-1 matrix    ( ⊺ scaled by a Dirac- distribution along the edge.The force balancing condition (1) translates to the statics equation while the length objective (2) relaxes to minimize where  1 (),  2 () are the eigenvalues (principal stresses) of , and 0 < ,  < 1.In the limit of ,  → 0, we have ( 2 =1 |  ()|  ) 1/ → rank(), and the integral approaches the size 2 of the region where rank() is nonzero, recovering (2).In the other limit where  =  = 1, (4) recovers the variational problem of [Strang and Kohn 1983] for  1 -minimal stress distributions.
In the rest of this section, we derive (4) from ( 2) and analyze the involved mathematical operators.We first describe a planar graph as a distribution over the plane, invoking the notion of a varifold.Next, we represent the varifold algebraically by a symmetric matrix field using the Veronese map.After this process, the notion of stress tensor field in (4) naturally emerges.Finally, the , -norms of (4) are justified using compressive sensing techniques.

𝜌
In the language of geometric measure theory, any planar weighted graph is a 1dimensional weighted varifold in the plane, which is a signed measure over the Grassman- which is the space of all unsigned planar directions.Note that G(1, R 2 ) ≃ S 1 .Intuitively, a measure over the space R 2 × G(1, R 2 ) of position and orientation indicates the occupancy of the graph.We denote the measure (i.e.weighted varifold) corresponding to our weighted graph . The total length of the edges in the planar embedding of this graph is given by the 1-dimensional volume of the support of this measure, i.e., (5) Here supp() is considered as an indicator function, and H  denotes the -dimensional Hausdorff measure.The problem can now be viewed as searching a measure over the nonlinear manifold R 2 × G(1, R 2 ).Next, we work with the varifold linear algebraically by invoking the Veronese map Sym .That is, V maps each orientation to a rank-1 symmetric matrix by the outer product.This Veronese map induces a pushforward map for measures.In particular, we have a linear map from weighted varifolds to symmetric-matrix-valued measures over the plane: e  e ⊺   x .Physically, V #  represents the stress distribution of the weighted varifold .On the space of symmetric matrix measures, we have the divergence operator div : given by contracting the differential and the tensor along one index: This is the same as the standard tensor divergence frequently used in elasticity.The divergence of the symmetric tensor represents the net traction force from the stress.In other words, the stress equilibrium condition is the varifold equivalent of (1).Here, ) represents the force distribution as a vector-valued measure.Using the above representation, we find (2) equivalent to the optimization problem minimize We recognize this as an instance of the sparse basis pursuit problem, which aims to find a sparse linear combination of atoms in a dictionary which satisfy underdetermined linear constraints [Chen et al. 2001].Here, the atom set consists of the point measures , and the objective on the support of our measure encourages sparsity with respect to this basis.
By performing a change of variables  = V #  via the linear map (6), Problem (9) is equivalent to finding a symmetric stress tensor field In compressive sensing, a standard approach converts such problems to norm minimization subject to linear constraints [Donoho and Elad 2003].Namely, define a linearly homogeneous function3 (referred to as a "norm") ∥ • ∥ whose unit ball has sharp corners containing the atom set (Fig. 3).Under linear-constrained norm minimization, these sharp corners generally lead to sparse solutions, as it is likely that the constraint affine plane will be tangential to the norm levelsets at these corners, which consist of few atoms.One can design such a unit ball by, for example, taking the convex hull of the atom set; this choice corresponds to the  1 convex relaxation and coincides with the minimal stress problem of [Strang and Kohn 1983].However, as we demonstrate in Section 5.1.5,this convex relaxation is too crude to produce sparse stress distribution in a general setting.Even sparser results can be obtained by using sharper star shapes, sacrificing convexity.Our method adopts the (nonconvex) (, )spectral norm with 0 <  < 1, 0 <  < 1.With the singular value decomposition  = ( Σ ⊺ ) , the spectral -norm4 is the -norm applied to the singular values Finally, our optimization problem takes the form (cf. ( 4)) minimize subject to div  = f . (11)

Divergence and the Killing Operator
Here, we analyze the divergence operator (7) for symmetric tensor measures, as it is the linear operator central to our optimization problem (11).The divergence operator is a linear operator from the space of symmetric tensor measures to the space of vector measures.Its adjoint is minus the Killing operator K [Berger and Ebin 1969;de Goes et al. 2014] completing the canonical duality diagram: The dual pairing between an element ) is given by 1 /2 of the integrated Frobenius pairing of their matrix representations and the dual pairing between f ∈ Γ( R 2 ⊗ ∧ 2  * R 2 ) and  ∈ Γ( * R 2 ) is given by Using this establishment of the duality relation, we can describe the solvability of the linear constraint in ( 11) The linear constraint ( 16) admits a solution  if and only if f ∈ im(div) = ker(K) ⊥ where (•) ⊥ denotes the annihilator space.Note that the kernel of the operator K is the collection of the one-forms corresponding to Killing vector fields (via ♭ R 2 ), which generate isometric flows on the domain.On R 2 , these Killing vector fields are the generators for rigid body transformations.In particular, ker(K) = span{, ,  −  }, where (, ) = ( 1 ,  2 ).Therefore, we have the following characterization for the valid f for ( 16).
Theorem 1. Eq. ( 16) admits a solution  if and only if the prescribed force distribution f =   e   satisfies the conditions of vanishing net force and vanishing total torque In our problem, we assume that the prescribed force distribution f satisfies the physically intuitive necessary and sufficient conditions (17).Even if we are handed a force distribution f that violates (17), it is straightforward to project it to fulfill the conditions by adding a suitable rigid motion vector field.

Representing Tensors and Differential Operators
We can represent symmetric tensors as arrays of their matrix elements.For instance, in the 2 × 2 case, a tensor with elements    can be represented as a vector ( 11 ,  22 ,  12 ) ⊺ ∈ R 3 .In Figure 4, we use this R 3 representation to depict the unit ball under our spectral -norm as well as the image of G(1, R 2 ) under the Veronese map V.By representing a symmetric matrix using these coordinates, the dual pairing (14) becomes and the differential operators div and K can be written as matrices of differential operators:

Obstacles
An extension to our optimization problem (11) allows us to construct varifolds which avoid "obstacles" placed in the domain.Let  : R 2 → R >0 be a function over the base space R 2 which takes the value 1 on obstacle-free areas of the domain, and takes the value  ≫ 1 on areas filled with obstacle.Then the objective minimize encourages solutions to avoid obstacles.Given large enough , solutions tend to pass around the obstacle-filled regions.

ALGORITHM
In this section, we describe an optimization algorithm equivalent to applying the Linearized Augmented Lagrangian Method [Yang and Yuan 2013] to our problem.This algorithm is best described as a Backward Euler discretization of a continuous gradient flow that optimizes the objective.The backward Euler steps are further translated into a sequence of variational problems using the method of incremental potential [Martin et al. 2011;Bouaziz et al. 2014;Li et al. 2020].
The aforementioned constrained optimization problem (20) takes the minimax form where the type of the Lagrange multiplier  is Γ( * R 2 ).The optimal solution can be found by following the coupled gradient descent and ascent flows with respect to  and .These continuous gradient flows are where  1 ,  2 are linear operators describing metrics for the space Γ( R 2 ⊗ sym  R 2 ⊗ ∧ 2  * R 2 ) of primal variables  and the space Γ( * R 2 ) of Lagrange multipliers , respectively.These metrics will be chosen suitably later during our derivation.
We discretize the flow (22) temporally using the backward Euler method.Replace the time derivatives with a step size , evaluate the right-hand sides of ( 22) at time step (•) (+1) , and approximate  (+1) on the right-hand side by  (+1) ≈  () +  () using  () =  () −  (−1) from the previous step to avoid a joint root finding system.The equations become: Let which simplifies (23a) into Solving ( 25) is equivalent to solving the optimization problem We choose  1 as the  2 Frobenius metric ∥ ∥ 2 Effectively,  1 is the identity map on the tensor coefficients; in particular, we can omit  1 in (23a) and ( 24).With this choice of  2 Frobenius metric, the sub-optimization problem ( 26) has an explicit solution given as a local shrinkage step on singular values [Yang and Yuan 2013]: Let  Σ ⊺ be the singular value decomposition for  (+1) ; then the solution to ( 26) is Finally, to complete ( 23) we choose the metric  2 .Observe that combining ( 24) and ( 23b) yields an expression that involves adding to  () by a term  2 K −1 2 div  () .Note that K and div are unbounded operators since they differential operators, preventing one from using larger time steps .To better precondition this update,  2 should be chosen to "cancel out" the derivatives K and div.A good choice is  2 = − div •K =  div div * for any scale factor  > 0. With this choice (K −1 2 div) becomes a bounded linear operator.
Spectral method for K and div.To invert and apply differential operators like div and K to tensors, we can take the Fast Fourier Transform (FFT) to these tensors defined on a rectangular domain of size  1 ×  2 with periodic boundary condition.FFT converts each partial derivative   to i 2    , where   is the integer index in the Fourier domain in the -th direction.All div and K operations (19) as well as (− div K) −1 can thus be performed in the Fourier domain as frequency-wise small complex matrix multiplications.Namely, applying the Fourier transform to the matrix representations in ( 19) gives the complex matrices Frequency-wise multiplication of these matrices gives a discrete representation of the operators div and K that can be used to implement Alg. 1. Explicitly, we follow the procedure listed in Alg. 2. An artifact of this procedure is the introduction of periodic boundary conditions on our domain; however, these can be nullified by placing boundary obstacles using our method.

RESULTS
We implement5 our algorithm in VEX in Houdini FX 18.5.All results were computed using a grid with resolution 256 × 256, a constant large time step  =  = 30, a fixed exponent  =  = 1 /2 (except for an ablation study), and 500 iterations.The computation were performed on a 2019 MacBook Air using a 1.6 GHz Dual-Core Intel Core i5.Each iteration takes about 40 ms with the FFT being the main bottleneck.6That is, each example is obtained within 20 seconds.

tension compression
The visualization of the results shows the trace of the stress tensor  at each grid cell.When  at a cell is rank 1 (i.e. it is in the image of the Veronese map), this trace shows the value of the only nonzero eigenvalue of .Blue strands represent tension, and can be interpreted as cables; red strands represent compression, interpreted as struts.

Numerical Tests, Validation, and Ablation
5.1.1Single cable.Consider a simple setup where the force distribution is given by two impulses on two points, and the force vector points radially away from each other (Fig. 5).The optimal support structure mediating this force distribution is the line segment joining the two points, representing a cable pulled by the forces.Fig. 5 shows the result of our algorithm over a few iterations, demonstrating that our method successfully reproduces a sharp line segment connecting the given points.
5.1.2Force sheets.Instead of prescribing forces concentrated on point sets, we test our algorithm for force distributed over onedimensional sets.Fig. 9 shows the result of our algorithm when the force is evenly distributed over two sheets facing each other.Despite Fig. 6.Sparse support structures found by our algorithm.Left: The force is set to evenly distribute over two sheets; this is the same result shown in Fig. 9 but the color axis is rescaled to display detailed support structures.Right: The same force distribution but with an additional disk obstacle between the force sheets; the algorithm avoids the obstacle and constructs "tied arch bridges" to maintain the stress-free condition.
the denseness of the force distribution, our method is able to find a sparse network of cables and struts supporting the given load.After a tone mapping, Fig. 6 (left) reveals the emergent detailed branches and reinforcement structure.We also test the algorithm for a more challenging obstacle configuration.In Fig. 6 (right), we place a large obstacle in between the force sheets of Fig. 9.The algorithm automatically finds a shockingly sophisticated system of cables and struts to wrap around the obstacle, held together by emergent "tied arch bridges" reaching beyond the convex hull of the support of the force.
5.1.4Convergence.Fig. 12 shows a typical convergence plot of the loss function for our algorithm.Note that there is no procedure in the optimization such as line-search that would shrink the optimization step size .With a fixed step size, the plot reflects an asymptotic stability of the flow (22).As expected from our loss function which encourages low rank tensors, Fig. 10 shows that the optimized stress is rank 1 or 0 except for a sparse set of branching points.
5.1.5(, ) dependency.Fig. 11 shows the force sheet setup for several combinations of parameters  and  in our (, )-spectral norm.Small  enforces global sparsity, while small  enforces local low-rank quality.In the case when  =  = 1, our objective reduces to  1 -nuclear norm minimization, and our solutions do not achieve the same sparsity.For all other experiments, we choose  =  = 1 /2.

Miscellaneous Examples
5.2.1 Bridge Designs.We test that our algorithm is able to generate reasonable results for common engineering problems.Fig. 13 shows our algorithm's output for a force distribution on a horizontal "road" applying uniform downward force, with two point support forces a distance below the road.A rectangular obstacle is placed a small distance below the road, to allow for vehicles or pedestrians to pass under the resulting structure.Our algorithm identifies an arch of struts with attached cables, similar to a real-world bridge design.
A similar setup is tested in Fig. 14, with a longer road, a thinner obstacle, and upward point forces a further distance below the road.Under these conditions, our algorithm generates a more organic bridge structure.This result is also presented in Fig. 2. 8 shows a manual bridge design that uses this result as blueprint.
Fig. 15 demonstrates another setup with two roads stacked vertically.Our algorithm generates a hybrid arch-suspension bridge with one tall arch.

Cantilever Beam.
The cantilever beam is a common test case in topology optimization research.In Fig. 16, we design a similar setup by placing a weight force at the end of a beam, and balancing forces on the other end representing points where the beam attaches to a wall.Our algorithm develops a series of curved cables and struts which support the weight at the end of the beam.These results appear similar to outputs from standard topology optimization routines.We develop an optimization algorithm which employs geometric-measure-theoretic techniques to find a sparse network of (potentially non-manifold) curves connecting a prescribed force distribution over a given domain.We have shown that our algorithm handles complex force distributions with obstacles in the domain using a series of Fast Fourier Transform and shrinkage steps.With the Fast Fourier Transform in our algorithm being the most costly step, we obtain each result less than half a minute.Moreover, the algorithm is able to produce strikingly nontrivial stress networks potentially useful in computer-aided designs.For example, Fig. 8 demonstrates a bridge model designed using the result of our algorithm.

CONCLUSION AND DISCUSSION
Though each solution obtained from our algorithm provides a physically plausible design, the method still demands rigorous theoretical guarantees.We have not proven convergence properties for our non-convex optimization problem.It is likely that concrete statements about convergence can be made by analyzing our Sobolev gradient flow.Moreover, we have not shown the mechanical stability of our results in addition to the equilibrium condition.
We also note that our (, )-spectral norm minimization formulation is only a non-convex continuous relaxation of the original combinatorial graph minimization problem.This can lead to suboptimal solutions in comparison to known minimal graphs (Fig. 17 The scope of this paper is limited to 2 dimensions with a Euclidean metric.Extending the theory to respect arbitrary Riemannian metrics requires replacing the derivatives in (7b) and (12b) by covariant derivatives, compromising the diagonalizability by Fourier Transform.Additionally, a 3-dimensional version of our algorithm has yet to be experimented thoroughly.Higher dimensions introduce visualization and storage challenges.These challenges can likely be resolved using techniques of [Liu et al. 2018] or [Palmer et al. 2022].
Finally, our treatment's use of low-rank stress tensors as varifolds opens new avenues in geometric measure theory.We hope that this physical interpretation leads to further developments in computational applications of geometric measure theory.
Fig. 2. A bridge with a hybrid suspension-tied-arch support structure designed by a simple geometric optimization algorithm derived from geometric measure theory.Given a user prescribed force distribution to support and obstacles to avoid (top left), the algorithm efficiently finds a sparse geometric measure (a varifold) representing a sparse stress distribution balancing the force (bottom right), performing at 25 iterations per second without GPU acceleration.The blue and red colors visualize the sign of the trace of the stress tensor, indicating tension and compression respectively.
Fig. 3.When the unit ball has sharp corners containing the atom set, constrained norm minimization leads to sparse solutions.

Fig. 4 .
Fig. 4. The | • |  unit ball (grey), and the image of the Veronese map (red) in the 3-dimensional space of symmetric 2 × 2 matrices.Solutions tend towards the sharp rims of this unit ball.

Fig. 7 .
Fig. 7. Our algorithm finds a sparse support structure for a pair of point forces with an obstacle (shaded region) blocking in between the points (cf.Fig. 1).5.1.3Obstacle.Fig. 7 shows the setup of Section 5.1.1 with an additional obstacle placed in between the points, similar to the illustration Fig. 1.The obstacle is prescribed through a large weight  = 1000 in the obstacle and  = 1 elsewhere.The algorithm finds a realistic structure avoiding the obstacle.We also test the algorithm for a more challenging obstacle configuration.In Fig.6(right), we place a large obstacle in between the force sheets of Fig.9.The algorithm automatically finds a shockingly sophisticated system of cables and struts to wrap around the obstacle, held together by emergent "tied arch bridges" reaching beyond the convex hull of the support of the force.

Fig. 8 .
Fig. 8.A bridge design inspired by the result of Fig. 14.

500 Fig. 9 .
Fig. 9.The result of running our algorithm for 500 iterations on f evenly distributed on two sheets, pointing horizontally away from the center line.See also Fig. 6, left.

Fig. 11 .Fig. 12 .
Fig. 11.The result of our algorithm with the setup of Section 5.1.2for varying ,  in the objective.

Fig. 13 .
Fig.13.A tied arch bridge discovered by our algorithm.f points downward on a horizontal "road", and points northeast and northwest at the lower left and lower right points respectively.An obstacle is placed in the bottom center, to allow boats or highways to pass under the bridge.

Fig. 14 .
Fig. 14.A hybrid suspension-tied arch bridge discovered by our algorithm.The horizontally distributed downward force represents a road, and upward forces at the two lower points represent supports.An obstacle is placed in the bottom center to allow boats or highways to pass under the bridge.

Fig. 15 .
Fig. 15.A hybrid arch-suspension bridge discovered by our algorithm.f represents the weight from two decks of roads, and normal force from two supports.There is a rectangular obstacle in the bottom center.

Fig. 16 .
Fig.16.The cantilever beam problem (left) is a common topology optimization test case.Our result (right) reproduces a pattern similar to the iconic Michell structure (left).We believe our algorithm can be used as a lightweight preprocessing step for more costly topology optimization methods.

Fig. 17 .
Fig. 17.Effect of our continuous relaxation to the original graph minimization problem.The result of our algorithm (right) for 4 point forces, pointing at 60 • from the horizontal at each point, is a sparse graph that deviates from the classical minimal Steiner tree (left).