Efficient Computation of Redundancy Matrices for Moderately Redundant Truss and Frame Structures

Large statically indeterminate truss and frame structures exhibit complex load-bearing behavior, and redundancy matrices are helpful for their analysis and design. Depending on the task, the full redundancy matrix or only its diagonal entries are required. The standard computation procedure has a high computational effort. Many structures fall in the category of moderately redundant, i.e., the ratio of the statical indeterminacy to the number of all load-carrying modes of all elements is less one half. This paper proposes a closed-form expression for redundancy contributions that is computationally efficient for moderately redundant systems. The expression is derived via a factorization of the redundancy matrix that is based on singular value decomposition. Several examples illustrate the behavior of the method for increasing size of systems and, where applicable, for increasing degree of statical indeterminacy.


Introduction
Information about the system-inherent properties of structures provides valuable insights into their behavior.For a usable spatially discrete truss or frame structure, the degree of statical indeterminacy is a widely used system-inherent property in elastostatics.It is equal to the difference between the number of unknown force and moment quantities and the number of linearly independent equilibrium equations.The statical indeterminacy is commonly considered as one aggregated integer number characterizing the number of redundant constraints within the entire structure.Unfortunately, statical indeterminacy does not describe the contribution of each member of the structure to load-carrying properties.Information about the spatial distribution of the degree of statical indeterminacy, as well as information about different load-carrying mechanisms, can be obtained from the redundancy matrix.This information shows the variety of load paths in the structure, which can be later used in the analysis and design of the structure.Particularly, the importance of individual members can be concluded from diagonal entries of the redundancy matrix, also called here distribution of statical indeterminacy.
A concept of redundancy in statically indeterminate truss and frame structures was proposed on the basis of the analogy between Gaussian adjustment calculus in geodesy and structural mechanics (Linkwitz 1961;Bahndorf 1991;Ströbel 1997).Recently, the idea of the redundancy matrix was extended to kinematically indeterminate structures (Tibert 2005;Zhou et al. 2015;Chen et al. 2018) and to continuum mechanical framework (Gade et al. 2021).A more extensive literature overview on the calculation and application of the redundancy matrix can be found in the work by von Scheven et al. (2021).The major applications of redundancy matrix in the design of structures aim primarily at the reliability and robustness of structures (Frangopol and Curley 1987;Pandey and Barai 1997;Kou et al. 2017;Spyridis and Strauss 2020) and the quantification of imperfection sensitivity (Eriksson and Tibert 2006;Ströbel and Singer 2008).More recent applications include adaptability, actuator placement, and optimized control in adaptive structures (Wagner et al. 2018;Geiger et al. 2020;Maierhofer and Menges 2019).A detailed discussion how the redundancy matrix can be used in the design and analysis of structures is beyond the scope of this paper and can be found in the references given above.
During the design process of structures, many different variants are usually studied.In this case, the redundancy distribution in the structure can be used as an indicator of the internal constraint or robustness of the different variants.These values can be either used as direct feedback to interactive design changes by an engineer or in an automatic optimization process (Bahndorf 1991).In both cases, a fast calculation of the redundancy matrix is essential.Recently, an algorithm for efficient updates of redundancy matrices in truss and frame structures is proposed in Krake et al. (2022).Therein, in distinction to the present contribution, generic algebraic formulations due to various modifications like adding, removing, and exchanging elements are employed.However, especially for larger structures, the calculation of the redundancy matrix itself is time-consuming as it involves the inverse of a matrix whose dimensionality is given by the number of degrees of freedom in the system.Therefore, this paper addresses the problem of efficient computation of the redundancy matrix or only its diagonal entries for large truss and frame structures, enabling rapid feedback to an engineer or faster optimization loops.In detail, we propose • a decomposition of the redundancy matrix based on a singular value decomposition (SVD), • application of this decomposition for efficient computation of the redundancy matrix or its diagonal entries, • a rough range where the proposed formula provides a significant speedup, • a set of benchmarks to investigate speedup.
In Section 2, all relevant aspects of matrix structural analysis and the definition of the redundancy matrix are provided.In Section 3, the SVD-based decomposition of the redundancy matrix is introduced and applied to derive efficient algorithms for the computation of the entire redundancy matrix and its diagonal entries only.The examples in Section 4 demonstrate the application of the proposed efficient computation formulations for different types of structures.
Here, the dependency of the speedup on the relative degree of statical indeterminacy and the problem size is discussed.The paper concludes with a summary and an outlook.

Background
In this section, a brief summary of matrix structural analysis is given, including the equilibrium, material, and redundancy matrices and their relevant properties.

Matrix structural analysis
Since the 1960s, static analysis for discrete models of spatial truss and frame structures purified the matrix notation (Argyris and Scharpf 1969;Przemieniecki 1968), which is shortly described for displacement-based version in the following.
Given is a discrete structural model consisting of  degrees of freedom,  n nodes, and  e elements.Elements are distinguished by nodal degrees of freedom and the number of loadcarrying modes  m .In the case of a pin-joined plane or spatial truss structure, each element possesses one load-carrying mode (axial tension/compression), i.e.,  m = 1.In the case of a frame structure, the number of load-carrying modes is  m = 6 for a spatial and  m = 3 for plane 2-node beam element.In general, a mixture of truss and beam elements is allowed in a structural model, i.e.,  m may vary between the elements.Therefore, the index  q is introduced that defines the number of all load-carrying modes of all elements in the model ( q =  m  e in case of a pure truss or beam model).
Compatibility equations relate generalized displacements d ∈ R  to the generalized elastic strains e el ∈ R  q and the generalized pre-strains e 0 ∈ R  q e el = Ad − e 0 . (1) Here, A ∈ R  q × denotes the compatibility matrix.Material equations relate the generalized elastic strains to the generalized stress resultants s ∈ R  q according to with C ∈ R  q × q being the material matrix.It is diagonal with positive entries for truss and beam elements, as described in von Scheven et al. ( 2021) and shown for a 2-node straight 3D beam in Appendix A. The equilibrium equations relate the generalized stress resultants to the external using the equilibrium matrix A ⊤ ∈ R × q .The irreducible form of the structural equations w.r.t the vector of generalized displacements satisfies these three sets of Equations ( 1) to (3).The square matrix is called the elastic stiffness matrix.It is symmetric by construction due to the symmetry of C following from its diagonal shape.
Several formal conditions on dimensions and indices of the model are assumed throughout the paper.First, the considered structures are statically indeterminate with a degree of statical indeterminacy  s =  q −rank(A ⊤ ).This condition ensures a nontrivial distribution of redundancy in the structures.Second, the structures are kinematically determinate or in other words, usable, i.e., rank(A) =  (Pellegrino and Calladine 1986;Pellegrino 1990;Pellegrino 1993).This condition ensures full rank of the elastic stiffness matrix K.

Redundancy distribution
Using the notations and assumptions of the previous subsection, the concept of the redundancy (Bahndorf 1991;Ströbel 1997;von Scheven et al. 2021) is described in the following.
The redundancy is a property independent of the external loads; thus, f = 0 is assumed.Elimination of the generalized displacement vector d from the compatibility Equation (1) using Equation (4) yields with the redundancy matrix It is a square matrix with a dimension equal to the number of all load-carrying modes  q .The main-diagonal entries   provide the distribution of the degree of statical indeterminacy of the structure (Bahndorf 1991;Ströbel 1997) between the load-carrying modes such that  q    =  s .The defining property of the redundancy matrix R in Equation ( 6) is that it maps the generalized pre-strains e 0 onto negative elastic strains, i.e., it reveals the incompatible, stress-inducing part of the action of e 0 .In the case of a statically determinate structure, R is a zero matrix and any e 0 leads to e el = 0.
The computational complexity for the calculation of the whole redundancy matrix is dominated by the inverse of the stiffness matrix and two matrix-matrix multiplications.The former operation has cost depending on the bandwidth of the stiffness matrix, but it is limited from above by the value for a dense matrix of the same size O ( 3 ).The latter operation has cost O ( •  2 q ).Since  is typically proportional to  q , the complexity scales as the third power of the problem size.
If only the distribution of the degree of statical indeterminacy in the structure is required, the calculation can be performed only for the main-diagonal entries of the redundancy matrix in Equation ( 7).Defining a  ∈ R 1× as the -th row of the compatibility matrix A, the following canonical formula can be used: This simple and straightforward modification of the definition of the redundancy matrix in Equation ( 7) makes use of simple linear algebra manipulations only to compute the diagonal entries.Therefore, the computational complexity is reduced to O (2 •  q ).Nevertheless, the computation of the inverse is probably the most expensive one.This operation is limited by O ( 3 ).Equation ( 8) serves as the reference for our proposed efficient computation of the diagonal of the redundancy matrix.
Premultypling Equation ( 6) with the material matrix yields The symmetric matrix CR is the so-called self-stress matrix, and Equation ( 9) is the map of the generalized pre-strains to the generalized stress resultants.
For further algebraic and spectral properties of the redundancy matrices in a discrete framework as well as redundancy functions in a continuous framework, we refer to the articles von Scheven et al. ( 2021) and Gade et al. (2021).

Method
In this section, a new formulation for the computation of redundancy matrices in large truss and frame structures is presented.According to the previous section, the costly computation of the redundancy matrix R, see Equation ( 7), or its diagonal entries, see Equation ( 8), is dominated by inverting the  ×  elastic stiffness matrix K.This step becomes more and more challenging with a growing number of degrees of freedom .From a mathematical point of view, however, the redundancy matrix has a rank equal to the degree of statical indeterminacy: This difference compared to the total number of load-carrying modes in the structure can be defined as the relative degree of statical indeterminacy  =   / q = ( q − )/ q , which satisfies 0 <  ⩽ 1.An analogous relative degree is known in adjustment of geodetic networks (Krarup 2006, p. 301).For small and moderate values of , the conventional computation of the redundancy matrix (involving the inverse of an  ×  matrix) may be inefficient.Our method builds upon this fact and is described in the following subsections.In the first subsection, a decomposition of the redundancy matrix is introduced.Then, the subsequent subsection makes use of this decomposition and presents our algorithms to efficiently compute either the entire redundancy matrix or the diagonal of the redundancy matrix.

Decomposition of the redundancy matrix
While Pellegrino (1993) used the SVD of the equilibrium matrix A ⊤ to relate it to physical properties, we propose to perform the SVD on the matrix product C 1 2 A to obtain new insights into the computation of redundancy matrices.The matrix C 1 2 A is also considered within a matrix force method by Robinson and Haggenmacher (1970).Let the SVD of C 1 2 A ∈ R  q × be given by where U = [u 1 , . . ., u  , u +1 , . . ., u  q ] ∈ R  q × q and V = [v 1 , . . ., v  ] ∈ R × are orthogonal matrices, and is a diagonal matrix with S = diag( 1 , . . .,   ) ∈ R × and  1 , . . .,   > 0. Mathematically, the first  column vectors u 1 , . . ., u  of the matrix U belong to the image of the matrix product C Tkachuk et al.

Efficient computation of redundancy matrices for moderately redundant structures
In contrast, the last  q −  column vectors u +1 , . . ., u  q belong to the left kernel of the matrix product C 1 2 A or, in other words, to the kernel of (C 1 2 A) ⊤ , i.e., for  =  + 1, . . .,  q holds: Now, we can integrate the above decomposition of C 1 2 A into the redundancy matrix formula.To do this, our first step is to rearrange the redundancy matrix equation in Equation ( 7) as follows: Since the matrices U and V are orthogonal, i.e., U ⊤ = U −1 and V ⊤ = V −1 , the inverse in the above equation can be calculated via the following formula: With this, we can rearrange the redundancy matrix equation via Equation ( 14) and Equation ( 15) further: The above equation reveals the mathematical matrix structure of the redundancy matrix.In fact, the conventional computation of the redundancy matrix translates into a computation of an orthogonal basis for the image of the matrix product C 1 2 A. Both approaches mainly reflect a computational complexity that is characterized by the number of degrees of freedom .However, our goal is to obtain a formula whose computational complexity is characterized by the difference  q − .To achieve this, we use the relation which reflects the connection between orthogonal projections (onto the image of C 1 2 A and kernel of (C 1 2 A) ⊤ ).In this way, we obtain our final result of the redundancy matrix using Equation ( 16) and Equation ( 17): The above final formula shows that the redundancy matrix can be computed via the  q − column vectors [u +1 , . . ., u  q ] and simple algebraic manipulations (e.g., multiplication with a diagonal Tkachuk et al.
Efficient computation of redundancy matrices for moderately redundant structures ⊲ orthogonal basis for the kernel of (C Algorithm 1 Efficient computation of the entire redundancy matrix in large truss or frame structures.matrix and outer products).In particular, this computation has a computation complexity that is characterized by the difference  q − .Thus, the smaller the difference  q −  is, the more efficient the formula gets.The only numerical issue is to compute the vectors [u +1 , . . ., u  q ] efficiently, i.e., in other words, to efficiently compute an orthogonal basis for the kernel of (C 1 2 A) ⊤ .This aspect is also discussed in the next subsection.
The column vectors u 1 , . . ., u  , u +1 , . . ., u  q characterize not only the kernel and image of (C 1 2 A) ⊤ but also the spectral components of the redundancy matrix.In fact, the columns of C − 1 2 [u 1 , . . ., u  ] are eigenvectors of R to the eigenvalue 0 and the columns of C − 1 2 [u +1 , . . ., u  q ] are the eigenvectors to the eigenvalue 1.To prove it, we use Equation ( 18) and consider the product of the redundancy matrix with an arbitrary column vector Another way of interpreting Equation ( 18) is to relate the formula to the self-stress matrix CR (see Equation ( 9)).In fact, Equation ( 18) can be rearranged into This formula shows that the matrix C 1 2 [u +1 , . . ., u  q ] decomposes the self-stress matrix CR.This decomposition is an intermediate result in Algorithm 1 (see Line 5).Therefore, our formula also covers the efficient computation of the self-stress matrix.
Analogously to the redundancy matrix, there is a relationship of the column vectors u 1 , . . ., u  , u +1 , . . ., u  q to the spectral components of the self-stress matrix.While the columns of C − 1 2 [u 1 , . . ., u  ] are eigenvectors of CR to the eigenvalue 0, the columns of C − 1 2 [u +1 , . . ., u  q ] solve the generalized eigenvalue problem of CR relative to C to the eigenvalue 1.To realize this, Equation ( 19) is multiplied with C (from the left), i.e.,

CR(C
which proves the eigenvalue statement mentioned above.

Efficient computation of the redundancy matrix
In the previous subsection, the formula in Equation ( 18) for the redundancy matrix is presented that favours an efficient computation if the difference  q −  is small.Based on this formula, we derive two different algorithmic approaches to compute either the entire redundancy matrix or only the diagonal of the redundancy matrix.Both scenarios use the column vectors [u +1 , . . ., u  q ], which is an orthogonal basis for the kernel of (C 1 2 A) ⊤ .This is the first step in both algorithms: Algorithms 1 and 2.
Numerically, the computation of an orthogonal basis for the kernel of (C 1 2 A) ⊤ can be performed via QR decomposition or another optimized function.Since our implementation was done with MATLAB, we propose the use of the built-in QR decomposition function via the command [Q, ∼, ∼] = qr(•) (with this notation, the QR decomposition makes use of permutation matrices to reduce fill-in).Even though the resulting orthogonal matrix Q = [u 1 , . . ., u  , u +1 , . . ., u  q ] also consists of the image of (C 1 2 A) and the kernel of (C 1 2 A) ⊤ , to our knowledge, it is the most efficient method to extract the desired column vectors [u +1 , . . ., u  q ].
In the following, we present the further steps for each of the two scenarios.We start with the efficient computation of the entire redundancy matrix.
Efficient computation of the entire redundancy matrix Algorithm 1 presents the efficient computation of the entire redundancy matrix.As explained previously, the computation of the column vectors [u +1 , . . ., u  q ] is the first step of the algorithm and is represented in Lines 2ś3.The next three lines implement the main formula that is described in Equation ( 18): while Line 4 and 6 simply perform matrix multiplications with diagonal matrices, Line 5 is a sum of outer product.
The computational complexity of the algorithm is mainly characterized by the computation of an orthogonal basis for the kernel of (C 1 2 A) ⊤ (as the other operations are very simple).Depending on the numerical method, the extraction of the  q −  vectors can be very efficient.This leads to a fast computation of the entire redundancy matrix for small and moderate values .
Efficient computation of the diagonal of the redundancy matrix Algorithm 2 presents the efficient computation of the diagonal entries of the redundancy matrix.As explained in the beginning of this subsection, the computation of the column vectors [u +1 , . . ., u  q ] is the first step of the algorithm and represented in Lines 2ś3.Since our goal is to only extract the diagonal entries of the redundancy matrix, we can rearrange Equation (18) as where e  is the -th unit vector and u , is the -th entry of the -th vector.The resulting formula is very simple and can be expressed in a more compact way, that is, the sum of each row of the matrix [u +1 , . . ., , where • is the Hadamard product (element-wise product).This procedure is the second step in Algorithm 2 and is described in Lines 4ś5.The computational complexity of the algorithm is dominated by the computation of an orthogonal basis for the kernel of (C 1 2 A) ⊤ (as the other operations are almost negligible).Depending on the numerical method, the extraction of the  q −  vectors can be very efficient.This leads to an extremely fast computation of the diagonal entries of the redundancy matrix for small and moderate values of .

Examples
This section demonstrates the application of the proposed efficient computation formulations for both computations of the distribution of the degree of statical indeterminacy and computation of the full redundancy matrix.With the first example, using a generic truss structure, the dependency of the efficiency of the proposed formulations on the relative degree of statical indeterminacy  will be demonstrated.The last two examples will apply the efficient computation Tkachuk et al.Efficient computation of redundancy matrices for moderately redundant structures ⊲ orthogonal basis for the kernel of (C Compute sum of each row of [u +1 , . . ., u  q ] • [u +1 , . . ., u  q ] ⊲ Equation ( 22) 6 end function Algorithm 2 Efficient computation of the distribution of the degree of statical indeterminacy (diagonal of the redundancy matrix) in large truss or frame structures.
to a more realistic truss structure and a three-dimensional frame system.All three examples are scalable to demonstrate the efficiency of the proposed algorithms for different problem sizes.For all examples, even a random distribution of the stiffnesses of the bars does not affect the increase in efficiency.All computations were done with MATLAB R2022a on a machine with a 2.80 GHz Intel Core i7-1165G7 processor and 32 GB RAM.An implementation of the algorithms and all examples is publicly available on DaRUS (Tkachuk et al. 2023).

Truss cylinder
The first example uses a cylindrical truss structure as it has been used for gas holders.The cylinder has a radius of 1 m and a height of 10 m.To examine structures with different numbers of elements and degrees of freedom, the number of segments in the axial and circumferential direction is defined by  and varied between 5 and 100.Furthermore, the degree of statical indeterminacy is varied by using different numbers of diagonal bracings.The relative degree of statical indeterminacy  will thus vary between 0.1 and 0.4.Figure 1 shows three structures with  = 8 and  = 0.1, 0.25, 0.4.Young's modulus  and cross-sectional area  are constant for all elements.The nodes in the lower level (marked in red) are fixed in all three directions.
Three series of comparisons with different values of  were examined for both the computation of the full redundancy matrix and of the distribution of the degree of statical indeterminacy.For each series,  was varied from 5 to 100 by steps of 5.But due to memory limitations on the reference machine, the series for the computation of the full redundancy matrix could not be completed up to  = 100.
In Figure 2, the results are shown for the computation of the full redundancy matrix.The left diagram shows the absolute computation time in seconds for the standard algorithms and the proposed efficient computation depending on the number of segments  for the truss cylinder with  = 0.1.It can be observed that the absolute gain in computation time grows with .The relative speedup for all three values of  is shown in the right diagram.For  = 0.1 (purple line), the speedup is in the range between 4 and 6 for most problem sizes.For the configurations with a higher relative degree of statical indeterminacy , only the speedup is shown in the right diagram.It can be observed that the speedup highly depends on the value of .For this truss system and  = 0.4, the speedup goes down to values between 1 and 2. This means, the proposed efficient algorithm is still faster for all values of , but the benefit is not as large as for systems with a lower relative degree of statical indeterminacy .As a second scenario, we consider now the computation of the distribution of the degree of statical indeterminacy, i.e., only of the diagonal entries of the redundancy matrix.The proposed Algorithm 2 is compared against the reference method described by Equation (8).In Figure 3 the absolute computation times are visualized for the cylinder with  = 0.1 (left) and the speedup for all three values of  (right).In general, we can observe a similar behavior as for the computation of the full redundancy matrix.The absolute gain in computation time grows with  and the speedup decreases with the relative degree of statical indeterminacy .However, the speedup is generally at a higher level and even reaches values up to 10 for  = 0.1.

Mero roof
The second example examines the calculation of the distribution of the degree of statical indeterminacy and the full redundancy matrix for a three-dimensional curved MERO space truss, as it is used in many buildings all over the world, e.g., in the Stockholm Globe Arena or the new Leipzig Trade Fair.The structure as shown in Figure 4 in two perspectives consists of two offset layers of square cells connected with diagonal bars.Each cell is 1 m × 1 m in the  −  plane.In this example, a quadratic function is used for the curvature in  and  direction.Young's modulus  and cross-sectional area  are constant for all elements.The nodes at the corners of the bottom layer (marked in red) are fixed in all three directions.For this structure the total degree of statical indeterminacy is equal to 45.The spatial distribution of the degree of statical indeterminacy in the structure is shown in Figure 5. Elements with a small contribution to the degree of statical indeterminacy are drawn with thin lines in light blue, while elements with a larger contribution are shown as thicker lines in darker blue.It can be observed that, e.g. the three bars connected to the top corners have a very small contribution to the degree of statical indeterminacy.In fact, they are statically determinate and therefore have a value of zero.In contrast, the bars in the lower layer generally have a higher contribution to the degree of statical indeterminacy as these bars undergo more constraint due to the four supports.To examine different problem sizes, the number of cells  in  and  direction is varied between 5 and 85.Since the topology of the structure is given, the relative degree of statical indeterminacy cannot be varied.For larger values of ,  tends toward the value 0.24.
In Figure 6, the computation time and speedup are shown for the full redundancy matrix (left) and the distribution of the degree of statical indeterminacy (right).Due to memory limitations, the full redundancy matrix could only be calculated up to  = 60.For both cases, it can be observed that the computation time for the proposed efficient computation is much smaller compared to the standard computation.The speedup for the computation of the full redundancy matrix is around 3 and for the distribution of the degree of statical indeterminacy even up to 5 or 6.
This example shows that the proposed efficient computation allows a significant acceleration of the computation even for real spatial trusses with more than 50,000 elements.

3D frame system
In the last example, efficient computation is applied to a three-dimensional frame structure.The structure shown in Figure 7 is a hyperbolic paraboloidal grid shell with square cells.Two neighboring edges (marked in red) are clamped to support the structure.All beams are rigidly connected and have the same Young's modulus  and cross-section.For spatial beam elements, the number of load-carrying modes is  m = 6.Therefore, also the degree of statical indeterminacy for every element can be up to 6.For this hyperbolic paraboloidal grid shell, the relative degree of statical indeterminacy tends toward  = 0.5 for large values of .In Figure 8, the computation time and the speedup are shown for the full redundancy matrix (left) and the computation of the distribution of the degree of statical indeterminacy (right).Also, for this more complex mechanical model and a higher value of the relative degree of statical indeterminacy, a considerable speedup can be observed.For larger values of , the speedup obtained by the proposed efficient computation for the distribution of the degree of statical indeterminacy is around 8 and for the efficient computation of the full redundancy matrix still between 3 and 4.

Conclusion
The paper addresses the problem of efficient computation of redundancy matrices for large truss and frame structures.Our proposed new formula is based on a SVD decomposition of the matrix product between the square root of the material matrix and the compatibility matrix.The few last columns of the decomposition construct the eigenspace of the redundancy matrix and provide a compact expression for the redundancy matrix.Our experiments have shown that for moderate values of a relative degree of statical indeterminacy, i.e., around  = 0.4, a speedup is observed for the full redundancy matrix computation.For lower values of , our proposed method allows for much higher speedups (up to 5 times faster).This speedup is stable over the number of members in a structure, as the three considered examples with scalable structures illustrate.Even a more significant speedup (up to 10 times faster) is observed for the computation of the distribution of the degree of statical indeterminacy, which is given by diagonal entries of the redundancy matrix.
The current study is limited to computations on CPU and straight structural members.In the future, hardware acceleration for the new computation algorithms on GPU may be considered.Curved beams are usually discretized with elements with more nodes and larger bandwidths.Thus, the behavior of the proposed algorithm may differ for frames with curved members.Finally, a new update algorithm for the full redundancy matrix based on our proposed decomposition can be developed and compared with the existing update algorithm from Krake et al. (2022).The implementation of the update should be straightforward, e.g., it profits from existing algorithms for dense and sparse matrices (Reichel and Gragg 1990;Andrew and Dingle 2014).
A Factorization of the stiffness matrix for a 3D Euler-Bernoulli beam element Herein, the expression from (Przemieniecki 1968) for the elastic stiffness matrix of a 2-node straight Euler-Bernoulli beam in 3D with a constant cross-section along the length is used.The local coordinate system { ê } 3 =1 is aligned with the axial direction for  = 1 and the principle bending axes  = 2, 3 with the second moments of inertia   and   , respectively.A factorization of the elastic stiffness matrix K  with a diagonal material matrix reads C  = 1  diag    3    3    , (A.1) where  is the cross-sectional area,  is the torsional moment of inertia,  and  are Young's and shear moduli, respectively.0 in the latter expression are zero column-vectors 0 3×1 .The load-carrying modes are axial tension/compression, torsion, and two bending/shear modes along each principle bending direction.As described in (von Scheven et al. 2021, pp. 45/46), the normalization is included in the matrix C  , so that K  = A ⊤  C  A  .This ensures a redundancy matrix that is independent of the chosen unit for the length dimension.

Figure 1
Figure1Truss cylinder: Systems with  = 8 and  = 0.1, 0.25, 0.4.Top view and perspective view.For a better visualization, the vertical axis is scaled differently than the axes in the horizontal plane.

Figure 2
Figure 2 Truss cylinder: Computation times for  = 0.1 and speedup for  = 0.1, 0.25, 0.4 for computation of full redundancy matrix.

Figure 3
Figure 3 Truss cylinder: Computation times for  = 0.1 and speedup for  = 0.1, 0.25, 0.4 for computation of distribution of the degree of statical indeterminacy.

Figure 5
Figure 5 Mero roof: Distribution of the degree of statical indeterminacy for  = 6.

Figure 6
Figure 6 Mero roof: Computation times and speedup for computation of full redundancy matrix (left) and distribution of the degree of statical indeterminacy (right).

Figure 8
Figure 8 3D frame system: Computation times and speedup for computation of full redundancy matrix (left) and distribution of the degree of statical indeterminacy (right).