Geometric Stress Functions, Continuous and Discontinuous

In his work on stress functions Maxwell noted that given a planar truss the internal force distribution may be described by a piecewise linear, $C^0$ continuous version of the Airy stress function. Later Williams and McRobie proposed that one can consider planar moment-bearing frames, where the stress function need not be even $C^0$ continuous. The two authors also proposed a discontinuous stress function for the analysis of space-frames, which however suffers from incompleteness. This paper provides a discontinuous stress function for $n$-dimensional space frames that is complete and minimal, along with its derivation from an $n$-dimensional continuous stress function.


Introduction
In his work on stress functions Maxwell [1] noted that given a planar truss the internal force distribution may be described by a piecewise linear, C 0 continuous (polyhedral) version of the Airy stress function [2].The bar forces correspond to the change of slope of the polyhedral Airy stress function, as a discrete version of the second derivatives.The same paper also deals with graphic representation of force distributions.It relates the polyhedral Airy stress function with the reciprocal diagram (Cremona force-plan [3] rotated with 90 degrees) of the planar truss.It also gives a scalar valued stress function for spatial problems the discrete analogue of which is compatible with force-diagrams of spatial trusses according to the representational idea of Rankine [4]; although Maxwell himself shows this stress function is not a complete description of static equilibrium.He then gives a vector valued stress function for spatial problems that is a superposition of three orthogonal copies of Airy's planar description.This vector valued stress function was later shown to be complete [5] for most engineering purposes.
The marriage of graphic representation with stress functions was also heavily emphasised in some recent works in graphic statics: Williams and McRobie [6] proposed that one can consider planar moment-bearing frames, where the stress function need not be even C 0 continuous.The three dimensional version of this discontinuous function, intended for space-frames was also introduced [7,8] relying on the scalar valued stress function of Maxwell; where the authors themselves have shown their stress function to be an incomplete description of static equilibrium.It was later shown that among single layer grid-shells what this description can handle are discretized structures of membrane shells that have Airy stress functions linearly proportional to their shape functions [9].As there are perfectly good grid-shells that don't satisfy this constraint we believe a different, less restrictive discontinuous stress function will be a useful addition to the literature.
This paper provides a dimension-independent discontinuous stress function for space-frames that is complete and minimal.On the application-side this allows for the efficient computation of spaceframes without the method prescribing spatial constraints on the structure.On the theoretical-side we believe the derivation of the stress function has additional explaining power and as such we present it here.During the derivation process we were guided by three requirements we set observing past works.
The first two requirements should be self explanatory.The third is motivated by the desire to solve the problem entirely, to give a description that does not loose its explaining power as the dimension of the space grows.When describing spatial rotations Weyl [10] phrased this rather eloquently: "The above treatment of the problem of rotation may, in contradistinction to the usual method, be transposed, word for word, from three-dimensional space to multi-dimensional spaces.This is, indeed, irrelevant in practice.On the other hand, the fact that we have freed ourselves from the limitation to a definite dimensional number and that we have formulated physical laws in such a way that the dimensional number appears accidental in them, gives us an assurance that we have succeeded fully in grasping them mathematically." The three requirements set us up for a long walk on the border of mathematics and mechanics.We will need a considerable amount of mathematical tools which we introduce below.

Notation, preliminaries
We will use k-vector valued differential forms.One way of looking at it is that one has the Grassmann algebra in R n where the scalar coordinates are replaced with differential forms.(Scalar valued forms in this interpretation, and the degree of them has to be the same in all coordinates).Real-number multiplication of the coordinates is replaced with the wedge product of differential forms.We will omit the wedge-product sign if possible and denote the "directions" of the Grassmann algebra with xi, xixj, . . .(for 1-vectors, 2-vectors, . . . ) while the components of the differential forms with di, didj and so on.The coordinates (functions) will be labelled by upper indices corresponding to the indexsets, for instance a 3-form coordinate of a 2-vector will have components like α 12,123 d1d2d3 x1x2.We will mostly follow the convention to list xi and di in lexicographic order.We will need the Hodge-duals [11] of both k-vectors and k-forms (we assume the usual Euclidean metric).We denote the Hodge-star on k-vectors with ⋆x, mapping a k-vector to an (n − k)-vector.This map may be given on the base vectors as ⋆x : where (1 . . .n) is the disjoint union of I and J and τ is the number of permutations required to bring (I, J ) to (1 . . .n). Similarly we have the one on scalar valued differential forms as where (1 . . .n) is the disjoint union of I and J and τ is the number of permutations required to bring (I, J ) to (1 . . .n).Under the Hodge-dual of a k-vector valued differential form we will mean the form achieved by applying composition of the two stars, i.e: This way α ∧ ⋆α gives an n-vector valued n-form.
We will need the scalar product of k-vector valued m-forms α and β, defined as the output of which is a real number.
Forces will correspond to 1-vectors while moments to 2-vectors.Force F = F i xi acting at point x = x i xi has moment x ∧ F with respect to the origin.In R n this is a 2-vector with n 2 coordinates, as there exists a moment with respect to the ortho-complement of each plane.(A moment is a rotating effect of a force and rotations can happen in each plane of the space.)The moment introduced this way differs from the engineering moment in R 3 having the sign of the second component flipped.This is captured in the relations Stresses in general are (n − 1)-forms, that need to be integrated along n − 1 dimensional hypersurfaces.As a consequence they will be measured in [N/mm d−1 ].As an example, we express the x1 directional stresses in R 3 as which corresponds to components of the Cauchy stress tensor with the usual x, y, z description as where the minus sign comes from the lexicographic ordering (1d1d3 = −1d3d1).
Strains may be represented with a 1-vector valued 1-form ϵ.Given some volume V of the material the work of the stresses on the strains inside V may be expressed by integrating the volume form which is an n-vector, thus isomorphic to a scalar.Here ⋆xϵ is an (n − 1)-vector valued 1-form.The spaces of 1-vectors and (n − 1) vectors play the role of the usual vector-space and the dual space of linear functionals.Since in case of R n they are isomorphic to each other, no attempt is made here to assign the role of primal-space and dual-space.Body forces (self-weight) correspond to vector-valued n-forms, we will denote them with ρ.
The exterior derivative is taken coordinate-wise on the k-vectors, we will denote it with d( ), omitting the parentheses if no confusion arises.The exterrior derivative satisfies d 2 α = d(dα) = 0 for all twice differentiable forms.It follows from the computational rules that given k-vector valued forms α and β the Leibniz Rule holds, where deg(α) denotes the degree of α.We will heavily rely on the following two results [12]: Lemma 1 (Poincaré Lemma).If dα = 0 throughout a simply connected region, then ∃β : α = dβ.

The continuous stress function
Let us cut out some volume V of the material.The equilibrium of the stresses acting on its boundary surface ∂V and of the body forces acting on V may be expressed as which must hold for all possible V , thus we have We will be able to give a stress function if there exists a potential function π for ρ such that ρ = dπ, as then d(σ + π) = 0 holds, implying the existence of (n − 2)-form ψ such that σ + π = dψ.Neither ψ nor π is unique, so much so that we may prescribe some of their components to be 0. We will treat ψ here and we will return to π later.It will be apparent later that we need where P is an index-set of n − 1 elements.To see that this is possible, assume we have found α such that σ + π = dα.We may create (n − 3)-form λ as where Q is an index-set of n − 2 elements, τ is the number of permutations required to bring the index-set (i, Q) to lexicographic order and dxi denotes integration with respect to the xi direction.
For n = 2 both P and Q is empty and any 0-form ψ is good.
For n = 3 only Q is empty.As an example the x2 direction of α looks like the undesirable part being α 2,2 .By integrating it we get function λ 2 = α 2,2 dx2.The stress and potential components pointing in the x2 direction will be calculated as which can be compared with Equations ( 14) and ( 15) are the same if Both are satisfied since we have α ∈ C 2 due to σ ∈ C 1 , for Equation ( 9) to make sense.We can see from the definition that the integration is always with respect to a single variable and is always possible without having to solve a system of differential equations.
The moment of the stresses acting on a small piece of hyper-surface at location x ∈ R n is calculated as x ∧ σ (with respect to the origin of the coordinate system).Similarly the moment of the body-forces may be expressed as x ∧ ρ.Cutting out some volume V of the material, the equilibrium of moments may be expressed as for any volume V .Since We already know d(σ + π) = 0, we have (here dx = i=1...n 1 dixi).At this point we have to go back to the fact that π is not unique and prescribe similarly to ψ.This guarantees dx ∧ π = 0 and we may conclude: The question becomes: Can we reconstruct ψ from dω and if so how?For an arbitrary differential form α the map α → dx ∧ α contains information loss, but for certain forms it is actually reversible.The idea can be seen in R 3 with the cross product as: Recalling how the cross product is the combination of the wedge product and the Hodge dual, the two maps become For differential forms ψ satisfying the orthogonality condition [ψ; dx] = 0 map (22) interchanges the indices as ψ P,i = (−1) n ψ i,P taking the vector valued (n − 2)-form to a (n − 2)-vector valued 1-form; furthermore map (23) is the inverse of map ( 22).The condition in Equation ( 10) is sufficient (but not necessary) to satisfy [ψ; dx] = 0, but getting rid of redundant parameters is useful in general, so let us parametrize the (n − 2)-form ω as follows: If α = dω then α ij,P = 0 ⇐ (i ∈ P and j ∈ P) and equivalently must hold.(For n = 2 ω is an arbitrary 0 − f orm.)In other words for any xixj moment component there is a single non-zero component ω ij,Q , exactly the one where (1 . . .n) is the disjoint union of (i, j) and Q.
To sum up the usage of what has been derived: Given ρ one has to find the potential function π satisfying Equation (18).Then choose any (n − 2)-form ω satisfying Equation (25) and the boundary conditions corresponding to the problem.The stresses are determined as (26)

Mechanical interpretation
We will be able to give a mechanical interpretation in case of ρ = 0 only, because in this case only does Equation ( 16) lead to d(x ∧ σ) = 0 implying the existence of ψ M , such that dψ M = x ∧ σ.
With this, we may consider a structure and cut in in half along a hyper-surface, denoting the surface of the cut by S. Then the moments of the stresses along S are given as We may note, that x where dδ = 0. Thus integrating the stress function ∂S ω provides a "correction term" to get correct moment values from the expression ∂S (x ∧ ψ).In case of n = 2 the stress function is actually a moment-function (a 0-form, having moment values over the plane).If the two endpoints of the cut are u, v ∈ R 2 , then ∂S ω = ω(v) − ω(u).(See for instance Phillips [13].)

Relation with earlier continuous stress functions
Using the notation of Sadd [14]: For n = 2 the stress function is the Airy stress function, ω = ϕ.
Here the Hodge dual of 2-vector x1x2 is a 0-vector, that acts like a scalar.The Hodge dual of a 2-vector valued 1-form is a scalar valued 1-form.

The discontinuous stress function
Rewriting the idea of Maxwell [1] in this language given a planar truss the internal force distribution may be described by a piecewise linear, C 0 continuous stress function ω, implying dψ = σ = 0 where there is no structure (between the rods) and dψ not being defined where there is structure.The planar case is somewhat degenerate as although dψ = 0 holds ψ cannot be the exterior derivative of anything since it is a 0-form; regardless, it is not hard to see that dψ = σ = 0 also implies the piecewise linearity of ω. (This treatment requires that the load of the structure is acting on joints as concentrated forces, ρ = 0 has to hold everywhere.)Williams and McRobie proposed [6] that one can consider planar moment-bearing frames, where the stress function need not be even C 0 continuous and have a well defined value at the axes of the rods.The method they proposed works in case of structures where the rod axes correspond to planar graphs, or equivalently the rod axes are edges of a planar mosaic.We will be able to actually strengthen this description to include structures with non planar-graphs, but for now we will assume the rod axes of the framework to correspond to edges (1-faces) of a convex, polyhedral n-dimensional mosaic.
If n ≥ 3, the non-existence of stresses outside the rod axes imply the existence of n − 3 form Ψ such that dΨ = ψ.This also gives us implying the existence of n − 3-form Ω such that γ = dΩ.Returning to the potential function in Equation ( 27), we may express the moment of the stresses as where we used d(x ∧ Ψ) = dx ∧ Ψ + x ∧ dΨ.This shape would only make sense if ψ and ω would be differentiable in the whole of ∂S.This will not be the case in general, we will have parts in which they are continuous, but the actual equilibrium will depend on what happens at the C 1 or even C 0 discontinuities.The use of this equation is that it tells us the shape of the stress function pieces on the continuous parts.Based on the above, when taking the discontinuous analogue of the continuous function we prescribe the following rules (not stricter than what has been done before) that will determine the discontinuous stress function: 1.The function need not be defined at the rod axes.
2. To every point outside the rod axes a single stress function piece has to correspond, satisfying dΨi = ψ i and dψ i = 0 on the entirety of R n .
We will start by looking at the n = 3 case below, then generalize.

The discontinuous stress function in 3 dimensions
Consider a tetrahedral piece of material, vertices of the tetrahedron denoted by p 1 , p 2 , p 3 , p 4 , point p being inside the tetrahedron (Figure 1).We want to replace this with 4 pieces of rods running from p i to p connected in a force and moment-bearing way.Choosing q 1 , q 2 , q 3 to be points close to p 1 on the edges of the tetrahedron, the force resultant of the stresses (in the continuous case) acting on the area enclosed by lines q i , q j may be calculated via line-integrals of 1-forms ψ ij on the respective lines as Here the three 1-forms are the same, we labelled them according to the curve-segments to introduce the logic of building up the resultant from parts.Distorting the geometry of the tetrahedron into the discontinuous case, we naturally get triangles p, p i , p j spanning 2-faces of the mosaic where the continuity of the line-integral may break.To express the condition σ = 0 we have to define the stress function (it has to be single-valued) in every point that is not on a 1-face (rod axis).We define a separate stress function piece corresponding to the inside of each 3-face and to the relative inside to every 2-face.(Relative inside of a k-face: points of the k-face that are not elements of a j-face such that j < k.)We don't really have a choice here, we cannot merge the pieces corresponding to the 2-faces since the 1-faces separate them.We cannot merge the pieces corresponding to the 3-faces with each other since the 2-faces separate them.If we merge the pieces of the 2-faces with the pieces of the 3-faces we get a single component that captures nothing of the structure.Thus the discrete version of Equation (32) will be where Ψij are the potential functions as introduced above, while ±Ψi are the discrete, orientation sensitive jumps corresponding to the line-integral passing through the 2-faces of the mosaic.(Here every unique label may denote a unique form.We avoid defining Ψi through the Dirac-delta function since we believe this to be a technicality.)Observing that the path-integral runs in a stress-free region since all the stresses are concentrated to the rod axis we can freely perturb the path without changing the resultant.Let us parametrize q i → p 1 with qi := (1 − t)q i + tp 1 where t ∈ [0, 1).We don't allow where the error term satisfies ϵi(t) → 0 as t → 1.Using this we may have where the error term ϵj(t) − ϵi(t) can be arbitrarily small.Thus we conclude Ψij(q j ) − Ψij(q i ) = 0 (36) must hold and Ψij cannot be used to parametrise the discontinuous case.We set these components to 0 on the basis that otherwise they would give us useless parameters.To properly treat the signs of ±Ψi(q i ) we assign an orientation to each 2-face, which can be captured by the normal vector ni (at the point of intersection with the path of integration).Let ti denote the tangent vector of the path of integration, at the point of intersection with the 2-face.The force in the bar (acting on the rod-star that the tetrahedron becomes, expressed in the global frame) will be where sign( ) denotes signum function and ⟨ni, ti⟩ is the usual scalar product.The requirement of sign consistency is not the distribution of ni but the fact that one chooses one and sticks to it through all the calculations.With this established, we may observe that in the discontinuous case we actually want Ψi to be constants, since we want the line-integral around the rod-axis to give the same force resultant even if the exact path is perturbed in the stress-free region.
We may similarly write up the moment of the force system with respect to the origin as The moment of the force system at p 1 with respect to p 1 is Thus to each 2-face corresponds a force system, and the stress function Ωi is −1 times the moment of the force system with respect to the points in the 2-face.As a consequence dΩi = dx ∧ Ψi (39) The integration is carried out on a tetrahedron containing point p 1 as an internal point.Components Ψ i correspond to the respective vertices q i of the tetrahedron while ψ is a surface integral to be evaluated everywhere except at the vertices.
holds (which is the discrete analogue of Equation ( 20)) and the force components stored in Ψi may be restored from the derivative, for instance in the way introduced in the continuous case.The sign sensitive sum can then be computed.Another way of looking at this is that any force system may be expressed as the vector pair (F i, M i) (the moment is again taken with respect to the origin.)The stress function is Ωi(x) = x ∧ F i − M i which makes sense on the whole of R 3 .The internal forces in the rod admit this decomposition as well.Thus Equations ( 37) and ( 38) can be expressed together as which is just a 6-dimensional vectorial sum.

Dimension-independent generalization
Given an n-dimensional mosaic with its polyhedral framework we have to decide to which elements of it do we attach a corresponding stress function piece.If n > 3 we actually have a choice in how to do this.Let us start by considering n = 4, where we can define pieces to each 2,3,4-face and try to merge them to minimize the components used.We cannot merge the components to the relative insides of 2-faces with each other (without involving higher dimensional faces) since the 1-faces separate them.We can however merge the pieces corresponding to the insides of 3-and 4-faces to a single piece.This is the case since in n = 4 the intersection of two 4-faces is a 3-face (or the empty set), and the internal points of the 3-face will be on the boundary of both 4-faces.In other words a 2 dimensional plane does not separate the 4 dimensional space, similarly how a line does not separate 3 dimensional space.With this we will have Ψi and Ωi corresponding to the two-faces as in case of n = 3 and ψ corresponding to the merged k-faces (k > 2).
We can have a similar set-up as in Figure 1, which lead to a line integral along a 2 dimensional simplex.In n = 4 we need to integrate along the boundary surface of a 3 dimensional simplex, that is a tetrahedron.Let us denote the tetrahedron with S, the point on the rod axis inside it with p 1 and the vertices of the tetrahedron with q 1 . . .q 4 (see Figure 2).The intersection of S with the two-faces of the mosaic are line segments p 1 q i .The resultant force acting on the rod may be expressed as ±Ψi(q i ). (41) Since we have d ψ = 0 we know ψ is finite in all points and a definite integral involving it does not change by removing a finite number of points from the domain.Thus and we again have a situation where the resultant is determined by the function components corresponding to the 2-faces and we set ψ = 0.This argument of merging pieces to get rid of anything but the pieces of the 2-faces works in any n > 3. Everything will work similarly to Equation ( 40), but we do have to generalize the sign convention.In general this may be done by choosing an orientation on each 2-face.If we pick point ci in it and place bi,1, bi,2 orthonormal base vectors there, the positive direction of rotation in the plane is represented by bi,1 ∧ bi,2 (see Figure 3).If we cut out an element of the rod, it will have two outwards pointing axial vectors a and −a.When determining the internal force system at the endpoints, the stress function will appear with positive sign if the respective axial vector ±a causes a positive directional rotation around ci.If we wish to express the internal forces of the structure at point r with outward normal a we may formally do this as where the bar is involved in m 2-faces.

Topological generalization, mechanical interpretation
At this point we can also consider what to do with non polyhedral frames.One effect of the orientation on the 2-face is that it can be identified with a direction of traverse of its boundary-loop.This loop exists if the 2-face is non planar, or the connectivity of the structure is more complicated.In the continuous case typically there is material in all points in the domain and there is static indeterminacy "everywhere".We define the continuous stress functions correspondingly "everywhere".For space frames the static indeterminacy is tied to the loops in the structure: one can not solve the static problem because the loop has no free end to start determining internal forces from.Thus we should define one stress function piece inside each "independent" loop, each loop introducing n + n 2 static indeterminacies.Thus we arrived at the observation that this is a combinatorical-linear algebra problem, where the combinatorical properties are determined by the topology of the structure.
These "independent" loops may be rigorously given by both algebraic topology and graph theory.The algebraic topological approach treats these loops as generators of the fundamental group of the graph of the structure, see for instance section 1.2 in the Algebraic Topology book of Hatcher [15].The graph theoretical approach treats these loops as generators of the cycle-space of the graph, see for instance section 1.9 in the Graph Theory book of Diestel [16].
The method given here can incorporate planar problems if they are embedded in at least R 3 , with potential functions Ω and Ψ having some constant 0 components.Furthermore, this description strengthens the existing planar one as non-planar graphs can also be computed this way.This will be seen in the proofs below.Before we present these proofs we take a look at the usage of what we derived.

Usage
The original idea of Maxwell turned loads and support reactions into internal forces by adding fictitious bars and at least one joint representing the "world" outside the structure, where the loads and support reactions meet.This way the equilibrium of the entire structure is expressed by the equilibrium of the added joint(s).We will refer to the larger structure created this way as the extended structure, to the non-extended one as the original structure.Since loads and support forces of the original structure will have to be sums of stress function pieces, this imposes conditions on the stress function pieces which can be expressed as a system of linear equations.In case the original structure is statically determinate there is a single solution to this system of linear equations.If the original structure is statically indeterminate we have to find the actual solution the structure chooses, depending on its geometry and material properties.If linear elastic materials and small displacements are assumed one possible solution strategy may be quadratic programming, to which we give a (3 dimensional) example below.

Example
Consider a tetrahedral frame, with moment bearing joints at points p 1 , p 2 , p 3 , p 4 !Let us load it with a concentrated force L ∈ R 3 at point p 1 !Support it with pinned supports at p 2 , p 3 and a roller at p 4 allowing only x3 directional force.The structure is drawn in the top of Figure 4. We will argue using the graph of the structure, joints at p i will correspond to graph vertices vi while a bar between p i and p j will correspond to edge {vi, vj}.To account for the loads and supports, we extend the graph by adding vertex v0 and additional graph edges {v0, vi} representing the load and the supporting forces.This is drawn in the lower part of Figure 4.
There are 10 loops in the graph, as follows: To each loop corresponds a stress function piece Ψi = (F i, M i), i = 1 . . .10, the coordinates of which will be our unknowns.We will adopt the notation that whenever referring to a force system present in a bar, we express the coordinates of the force system that acts on the joint with the smaller index.A stress function piece contributing to a bar force will appear with positive sign if the corresponding loop traverses the vertices of the edge corresponding to the bar in ascending order.As the load corresponds to edge {v0, v1}, the prescribed load turns into a condition on the stress function pieces as The condition that there are pinned supports at p 2 and p 3 mean the moment and force components satisfy Finally, the roller support at p 4 means where the last two equations represent directional constraint of the roller.Equations ( 44) -(49) may be collected as system of linear equations in the shape of where y contains the unknowns F i and M i, C is the coefficient matrix and b contains the effect of the load.The elastic deformational energy stored in the bar between p j and p k may expressed from the (F i, M i)-shaped dynames with the help of a 6-by 6 matrix [17], which may also be used to express this energy as a function of stress function coordinates.Rearranging these equations we may write up the total elastic deformational energy in the form E = j,k E j,k = 1 2 y T Qy (k > j) where Q = Q T and arrive at a quadratic programming problem in the form minmize:

Formal proofs
Unsurprisingly, since we derived the continuous stress function from the equilibrium equations and the discontinuous stress function from the continuous one, the discontinuous stress function is equivalent with the static equilibrium of the extended structure.We will show this in two steps.First we show that each internal stress distribution the stress function gives is automatically in equilibrium, then we show that any solution of static equations can be represented this way by a stress function that is unique (up to the choice of the global coordinate system).

Automatic equilibrium
Theorem 2. The proposed discontinuous stress function gives internal force systems that are in static equilibrium.
Proof.Consider a joint of the structure where j rods meet, and consider the corresponding vertex of the graph of the structure.Each loop of the graph that travels through the vertex enters on one graph-edge and exits on another.Thus when summing up the force resultants at the ends of the rods acting on the joint, each stress function component will be present twice, with opposite signs.For the equilibrium of forces we have (after some rearrangement) where k denotes the number of loops passing through the vertex.Similarly, we can write up the sum of the moments with respect to p, which after some rearrangement will take the shape of with the same indices, completing the proof.

Completeness and minimality
We will argue using the extended structure.(This is in contrast to the completeness investigation of continuous stress functions [1,18], where the question of completeness can not be investigated without considering what the boundary of the solid is like [19,5].Here we don't prescribe boundary conditions as they would only exclude certain loads and we are interested in parametrizing the general case.) Theorem 3. The proposed discontinuous stress function is a complete and minimal parametrization of the internal forces of the extended structure.
Proof.We have to show that whatever internal force distribution that is in equilibrium is given in the extended structure, there is a unique corresponding stress function built from the appropriate components.Recall how the force components can be calculated through component-wise summation (Equation (43)), where the topology of the structure determines the summations.We will manipulate a graph that will act as a topological aid to write up the correct equations.The starting shape of this graph will be the graph of the extended structure.We will calculate the stress function coordinates component-wise, row-by row.Let gi denote the x1 directional force component of each stress function piece, where i runs on all the generator-loops of the cycle space of the graph of the extended structure.Let f k denote the x1 directional component of the bar-force in bar k (k runs on all the bars).For each i in ascending order we may do the following: Find the loop in the graph corresponding to Ψi (Figure 5, left).Choose any bar k in the loop and express f k as where Ji is some index set.After this equation is written up contract the loop in the graph, unifying the involved vertices to a single one.By contracting the loop we make sure not to use component gi again.This way the index set Ji will satisfy: j ∈ Ji =⇒ j > i.
After we do this for all i we get a system of linear equations Ag = f , whose coefficient matrix A is square, upper triangular and each element on the main diagonal is ±1.The determinant of this matrix is ±1 (the product of the elements on the main diagonal), thus it is invertible and we may solve for g.We may also repeat the whole procedure for the other components of (F i, M i), in total n + n 2 times.As such we may find a suitable stress function distribution to any force system in the structure, that is in equilibrium.We may also note that the number of loops equals n + n 2 the degree of static indeterminacy of the extended structure, implying that in the general case of a frame we can not get away with less stress function parameters.
We still have to see, that choosing different bars in each loop does not give a different stressdistribution, or in other words f contains force components that parametrize the self-stresses of the structure.This can be seen by doing the loop-contraction procedure backwards.At each backwardsstep f k may be used as a parameter and the rest of the unknowns in the step may be calculated from  the static equilibrium equations (see Figure 6).At each backwards-step the number of added nodes equals the number of unknown components and there is an independent equilibrium equation corresponding to each joint of the structure (corresponding to each vertex of the graph).The mechanical interpretation of the equilibrium equation of a fictitious vertex (one that contains a loop contracted into it) is the sum of all the equilibrium equations of the joints that were present in the loop.As the procedure restores the internal force distribution of the structure, the proof is concluded.

Conclusion
Motivated by the previously open problem of finding a complete three dimensional discontinuous stress function we investigated the subject of stress functions in a systematic way.We based our approach on the differential-form nature of stresses, one of the cornerstones of the connection between elasticity and geometry.We rewrote the static equilibrium equations into a differential-form shaped dimension-independent continuous stress function, that in simply connected domains is equivalent with static equilibrium.Then, we took the defining property of idealized space-frames (stresses are zero everywhere except at rod axes) and applied it to our continuous function, thereby deriving a dimension-independent discontinuous stress function that is equivalent with static equilibrium of space-frames.This approach allowed us not only to solve the previously open problem, but we also improved on the planar construction of Maxwell by being able to treat planar mechanical problems with non planar graphs.Apart from this efficiency we could also see how and why the stress function description works: We saw that the stress functions are a parametrization of the self stresses of the structure and they should correspond to whatever is causing the static indeterminacy.In the continuous case the base-problem of elasticity is statically indeterminate, in the discontinuous case the roots of the in-determinacy are the loops in the extended structure.These loops are the generating elements of the first fundamental group of the extended structure, showing that the number of function-pieces required in the discontinuous stress function is determined by the topology and not the metric properties of the structure.Sticking with the idealized line-model of the structure and not taking material properties into account, these metric properties become important if one prescribes constraints in the internal force distribution, like introducing ball-joints enforcing truss-like behaviour.This will tie the discrete stress functions to line-geometry, as for special cases the dynames in Equation (43) turn into projective line-coordinates.We hope to continue this work by investigating space-trusses this way.We would not mind arriving at some graphic representation of the internal force distribution of space-trusses if possible, but we do wish to derive it from geometry instead of relying on a representation scheme rooted only in tradition.
Furthermore, although variational methods at first might seem far from geometry, the use of the discontinuous geometric stress function being equivalent with static equilibrium can be seen when using the Principle of the Minimum of Complementary Potential Energy.This principle requires one to take variations enforcing static equilibrium, which may be cumbersome if tried from a direct description of internal forces.Using the discontinuous stress function provided here one only has to take variations in the space of the stress functions, making the use of this principle trivial.

Remark 1 .
Admittedly, this is far from standard notation.The point of choosing it is to keep track of what object has what mechanical meaning.See for instance forces and moments below.

Figure 1 :
Figure 1: Tetrahedral piece of material and the discontinuous version of a line integral on it.The negative signs are due to the orientation of the boundary (endpoints) of the curve-segments.

Figure 2 :
Figure2: Decomposition of the surface integral in case of n = 4.The integration is carried out on a tetrahedron containing point p 1 as an internal point.Components Ψ i correspond to the respective vertices q i of the tetrahedron while ψ is a surface integral to be evaluated everywhere except at the vertices.

Figure 3 :
Figure3: Orientation of a 2-face of a mosaic may be expressed via the exterior product of two vectors in the mosaic-face.This implies a direction of traverse on it's boundary-loop.

Figure 4 :
Figure 4: An example to handle loads and supports as fictitious internal bars.The world outside the structure is represented by vertex v 0 , the equilibrium of the loads and supports correspond to the equilibrium of vertex v 0 .
Qy under constraint: Cy = b by relying on the Principle of the Minimum of Complementary Potential Energy.

Figure 5 :
Figure 5: Contracting loop i corresponding to equation i in the process of determining the stress function corresponding to a given force-distribution.

Figure 6 :
Figure 6: Doing the loop-contraction procedure backwards to determine non-parameter rod forces.