tests/test.sh exit code (0 → resolved); the classification below is post-hoc and cannot change it./opt/truss2d/src/solver.cpp
1 #include "truss2d/solver.hpp"
2
3 #include <cstddef>
4 #include <stdexcept>
5 #include <vector>
6
7 #include "truss2d/dense_matrix.hpp"
8 #include "truss2d/types.hpp"
9
10 namespace truss2d {
11
12 // ============================================================================
13 // TODO(candidate): Implement the static solver for the cable/strut/bar network
14 // described in instruction.md.
15 //
16 // The four member functions below are stubs that currently return zero-valued
17 // results, so the analysis is physically wrong and the test suite FAILS.
18 // Replace the stub bodies with correct implementations. Do NOT change the
19 // public signatures declared in include/truss2d/solver.hpp, and reuse the
20 // existing support layer (DenseMatrix, solve_spd, and the Model helpers
21 // element_length / element_direction).
22 //
23 // Read instruction.md for the governing model, output conventions, and
24 // sign/units conventions. This file is the ONLY file you need to edit.
25 // ============================================================================
26
27 DenseMatrix StaticSolver::element_stiffness_global(std::size_t e) const {
28 // TODO(candidate): single-member elastic 4x4 stiffness in the global frame.
29 (void)e;
30 return DenseMatrix(4, 4); // all zeros -> wrong
31 }
32
33 DenseMatrix StaticSolver::assemble_system(const std::vector<bool>& active) const {
34 // TODO(candidate): assemble active-member elastic stiffness + spring terms.
35 (void)active;
36 const std::size_t n = model_.num_dofs();
37 return DenseMatrix(n, n); // all zeros -> wrong
38 }
39
40 std::vector<double>
41 StaticSolver::prestrain_force(const std::vector<bool>& active) const {
42 // TODO(candidate): nodal force from member prestrain for this active set.
43 (void)active;
44 const std::size_t n = model_.num_dofs();
45 return std::vector<double>(n, 0.0); // all zeros -> wrong
46 }
47
48 SolveResult StaticSolver::solve() const {
49 // TODO(candidate): validate, resolve the active set, solve, recover
50 // displacements / reactions / axial forces / active flags.
51 const std::size_t n = model_.num_dofs();
52 const std::size_t ne = model_.num_elements();
53 SolveResult result;
54 result.displacements.assign(n, 0.0);
55 result.reactions.assign(n, 0.0);
56 result.axial_forces.assign(ne, 0.0);
57 result.utilization.assign(ne, 0.0);
58 result.active.assign(ne, true);
59 return result; // trivially zero -> wrong
60 }
61
62 } // namespace truss2d
63
/opt/truss2d/include/truss2d/solver.hpp
1 #ifndef TRUSS2D_SOLVER_HPP
2 #define TRUSS2D_SOLVER_HPP
3
4 #include <cstddef>
5 #include <vector>
6
7 #include "truss2d/dense_matrix.hpp"
8 #include "truss2d/model.hpp"
9
10 namespace truss2d {
11
12 /// Result of a static analysis of a cable/strut/bar network.
13 struct SolveResult {
14 /// Nodal displacements, length == model.num_dofs(), ordered by global DOF
15 /// (2*node + component).
16 std::vector<double> displacements;
17
18 /// Support reaction forces at every global DOF, length == model.num_dofs().
19 /// Non-zero only at DOFs of nodes carrying a spring support.
20 std::vector<double> reactions;
21
22 /// Internal axial force per element, length == model.num_elements().
23 /// Sign convention: positive = tension, negative = compression. A member
24 /// that has dropped out of the load path (a slack cable / a separated
25 /// strut) reports exactly 0.
26 std::vector<double> axial_forces;
27
28 /// Signed utilization per element, length == model.num_elements(). Inactive
29 /// unilateral members report 0.
30 std::vector<double> utilization;
31
32 /// Whether each member participates in the converged load path, length ==
33 /// model.num_elements(). A Bar is always active; a Cable is inactive when
34 /// slack; a Strut is inactive when separated.
35 std::vector<bool> active;
36 };
37
38 /// Static solver for a 2D network of two-force members (bars, cables, struts)
39 /// on spring supports, with member prestrain.
40 ///
41 /// The headline departures from a classical linear pin-jointed truss are
42 /// documented in the analysis contract (instruction.md); this header only
43 /// fixes the public surface. The helpers below expose the ordinary elastic
44 /// pieces used by callers and tests.
45 class StaticSolver {
46 public:
47 explicit StaticSolver(const Model& model) : model_(model) {}
48
49 /// Run the full analysis. Calls Model::validate() first and propagates any
50 /// exception it throws. Throws std::runtime_error if the governing system
51 /// is singular (e.g. a mechanism / under-supported network).
52 SolveResult solve() const;
53
54 /// The 4x4 single-member elastic stiffness contribution for element `e`
55 /// expressed in the global frame, local DOF order {n1.x, n1.y, n2.x, n2.y}.
56 /// This is the elastic part only (it does NOT encode kind, prestrain, or
57 /// supports). Exposed for testing/reuse.
58 DenseMatrix element_stiffness_global(std::size_t e) const;
59
60 /// Assemble the full governing stiffness matrix (num_dofs x num_dofs) for a
61 /// GIVEN set of participating members `active` (length num_elements()):
62 /// the sum of the elastic stiffness of every active member plus every
63 /// spring-support contribution. Inactive members contribute nothing.
64 /// Exposed so the equilibrium residual can be checked directly.
65 DenseMatrix assemble_system(const std::vector<bool>& active) const;
66
67 /// The global nodal force vector (length num_dofs()) produced by member
68 /// prestrain for a GIVEN active set, i.e. the prestrain contribution to the
69 /// right-hand side of the governing system. Inactive members contribute
70 /// nothing. Exposed so the equilibrium residual can be checked directly.
71 std::vector<double> prestrain_force(const std::vector<bool>& active) const;
72
73 private:
74 const Model& model_;
75 };
76
77 } // namespace truss2d
78
79 #endif // TRUSS2D_SOLVER_HPP
80
/opt/truss2d/include/truss2d/model.hpp
1 #ifndef TRUSS2D_MODEL_HPP
2 #define TRUSS2D_MODEL_HPP
3
4 #include <cstddef>
5 #include <vector>
6
7 #include "truss2d/types.hpp"
8
9 namespace truss2d {
10
11 /// A 2D structural network of two-force members with spring (skew-roller)
12 /// supports and per-member tension/compression character and prestrain.
13 ///
14 /// The model is a plain data container plus light validation / derived-quantity
15 /// helpers. It is fully implemented; the numerical analysis lives in
16 /// StaticSolver (see solver.hpp).
17 class Model {
18 public:
19 std::size_t add_node(double x, double y);
20
21 /// Add a member between two existing nodes.
22 /// `kind` - Bar (default), Cable, or Strut (see types.hpp).
23 /// `prestrain` - mechanical installed axial strain (default 0).
24 /// `alpha,dT` - optional thermal expansion data (defaults 0).
25 /// Throws std::out_of_range if a node index is invalid.
26 std::size_t add_element(std::size_t n1, std::size_t n2, Section section,
27 MemberKind kind = MemberKind::Bar,
28 double prestrain = 0.0,
29 double alpha = 0.0,
30 double dT = 0.0);
31
32 /// Add a grounded spring support at `node` resisting displacement along
33 /// direction (dx, dy) with stiffness `kappa` [N/m]. The direction is
34 /// normalized internally. Throws std::out_of_range if `node` is invalid,
35 /// std::runtime_error if (dx, dy) is the zero vector or kappa <= 0.
36 /// `settlement` (optional, default 0) prescribes the spring's grounded-end
37 /// offset along (dx, dy); the spring force depends on the node displacement
38 /// along that direction RELATIVE to this offset.
39 void add_spring(std::size_t node, double dx, double dy, double kappa,
40 double settlement = 0.0);
41
42 /// Apply a force `value` [N] to a single global DOF (2*node + component).
43 /// Throws std::out_of_range if the DOF is invalid.
44 void add_load(std::size_t dof, double value);
45
46 const std::vector<Node>& nodes() const { return nodes_; }
47 const std::vector<Element>& elements() const { return elements_; }
48 const std::vector<SpringSupport>& springs() const { return springs_; }
49 const std::vector<Load>& loads() const { return loads_; }
50
51 std::size_t num_nodes() const { return nodes_.size(); }
52 std::size_t num_elements() const { return elements_.size(); }
53 std::size_t num_dofs() const { return 2 * nodes_.size(); }
54
55 /// Undeformed length of element `e` [m].
56 /// Throws std::out_of_range if `e` is invalid.
57 double element_length(std::size_t e) const;
58
59 /// Direction cosines (cos, sin) of element `e` measured from node n1
60 /// toward node n2, in the global frame.
61 /// Throws std::out_of_range if `e` is invalid, std::runtime_error if the
62 /// element has zero length.
63 Vec2 element_direction(std::size_t e) const;
64
65 /// Throws std::runtime_error if the model is structurally ill-formed
66 /// (no nodes, no elements, a zero-length element, a member with
67 /// non-positive E or A, or no spring supports at all).
68 void validate() const;
69
70 private:
71 std::vector<Node> nodes_;
72 std::vector<Element> elements_;
73 std::vector<SpringSupport> springs_;
74 std::vector<Load> loads_;
75 };
76
77 } // namespace truss2d
78
79 #endif // TRUSS2D_MODEL_HPP
80
/opt/truss2d/include/truss2d/solver.hpp
1 #ifndef TRUSS2D_SOLVER_HPP
2 #define TRUSS2D_SOLVER_HPP
3
4 #include <cstddef>
5 #include <vector>
6
7 #include "truss2d/dense_matrix.hpp"
8 #include "truss2d/model.hpp"
9
10 namespace truss2d {
11
12 /// Result of a static analysis of a cable/strut/bar network.
13 struct SolveResult {
14 /// Nodal displacements, length == model.num_dofs(), ordered by global DOF
15 /// (2*node + component).
16 std::vector<double> displacements;
17
18 /// Support reaction forces at every global DOF, length == model.num_dofs().
19 /// Non-zero only at DOFs of nodes carrying a spring support.
20 std::vector<double> reactions;
21
22 /// Internal axial force per element, length == model.num_elements().
23 /// Sign convention: positive = tension, negative = compression. A member
24 /// that has dropped out of the load path (a slack cable / a separated
25 /// strut) reports exactly 0.
26 std::vector<double> axial_forces;
27
28 /// Signed utilization per element, length == model.num_elements(). Inactive
29 /// unilateral members report 0.
30 std::vector<double> utilization;
31
32 /// Whether each member participates in the converged load path, length ==
33 /// model.num_elements(). A Bar is always active; a Cable is inactive when
34 /// slack; a Strut is inactive when separated.
35 std::vector<bool> active;
36 };
37
38 /// Static solver for a 2D network of two-force members (bars, cables, struts)
39 /// on spring supports, with member prestrain.
40 ///
41 /// The headline departures from a classical linear pin-jointed truss are
42 /// documented in the analysis contract (instruction.md); this header only
43 /// fixes the public surface. The helpers below expose the ordinary elastic
44 /// pieces used by callers and tests.
45 class StaticSolver {
46 public:
47 explicit StaticSolver(const Model& model) : model_(model) {}
48
49 /// Run the full analysis. Calls Model::validate() first and propagates any
50 /// exception it throws. Throws std::runtime_error if the governing system
51 /// is singular (e.g. a mechanism / under-supported network).
52 SolveResult solve() const;
53
54 /// The 4x4 single-member elastic stiffness contribution for element `e`
55 /// expressed in the global frame, local DOF order {n1.x, n1.y, n2.x, n2.y}.
56 /// This is the elastic part only (it does NOT encode kind, prestrain, or
57 /// supports). Exposed for testing/reuse.
58 DenseMatrix element_stiffness_global(std::size_t e) const;
59
60 /// Assemble the full governing stiffness matrix (num_dofs x num_dofs) for a
61 /// GIVEN set of participating members `active` (length num_elements()):
62 /// the sum of the elastic stiffness of every active member plus every
63 /// spring-support contribution. Inactive members contribute nothing.
64 /// Exposed so the equilibrium residual can be checked directly.
65 DenseMatrix assemble_system(const std::vector<bool>& active) const;
66
67 /// The global nodal force vector (length num_dofs()) produced by member
68 /// prestrain for a GIVEN active set, i.e. the prestrain contribution to the
69 /// right-hand side of the governing system. Inactive members contribute
70 /// nothing. Exposed so the equilibrium residual can be checked directly.
71 std::vector<double> prestrain_force(const std::vector<bool>& active) const;
72
73 private:
74 const Model& model_;
75 };
76
77 } // namespace truss2d
78
79 #endif // TRUSS2D_SOLVER_HPP
80
/opt/truss2d/src/solver.cpp
1 #include "truss2d/solver.hpp"
2
3 #include <cstddef>
4 #include <stdexcept>
5 #include <vector>
6
7 #include "truss2d/dense_matrix.hpp"
8 #include "truss2d/types.hpp"
9
10 namespace truss2d {
11
12 // ============================================================================
13 // TODO(candidate): Implement the static solver for the cable/strut/bar network
14 // described in instruction.md.
15 //
16 // The four member functions below are stubs that currently return zero-valued
17 // results, so the analysis is physically wrong and the test suite FAILS.
18 // Replace the stub bodies with correct implementations. Do NOT change the
19 // public signatures declared in include/truss2d/solver.hpp, and reuse the
20 // existing support layer (DenseMatrix, solve_spd, and the Model helpers
21 // element_length / element_direction).
22 //
23 // Read instruction.md for the governing model, output conventions, and
24 // sign/units conventions. This file is the ONLY file you need to edit.
25 // ============================================================================
26
27 DenseMatrix StaticSolver::element_stiffness_global(std::size_t e) const {
28 // TODO(candidate): single-member elastic 4x4 stiffness in the global frame.
29 (void)e;
30 return DenseMatrix(4, 4); // all zeros -> wrong
31 }
32
33 DenseMatrix StaticSolver::assemble_system(const std::vector<bool>& active) const {
34 // TODO(candidate): assemble active-member elastic stiffness + spring terms.
35 (void)active;
36 const std::size_t n = model_.num_dofs();
37 return DenseMatrix(n, n); // all zeros -> wrong
38 }
39
40 std::vector<double>
41 StaticSolver::prestrain_force(const std::vector<bool>& active) const {
42 // TODO(candidate): nodal force from member prestrain for this active set.
43 (void)active;
44 const std::size_t n = model_.num_dofs();
45 return std::vector<double>(n, 0.0); // all zeros -> wrong
46 }
47
48 SolveResult StaticSolver::solve() const {
49 // TODO(candidate): validate, resolve the active set, solve, recover
50 // displacements / reactions / axial forces / active flags.
51 const std::size_t n = model_.num_dofs();
52 const std::size_t ne = model_.num_elements();
53 SolveResult result;
54 result.displacements.assign(n, 0.0);
55 result.reactions.assign(n, 0.0);
56 result.axial_forces.assign(ne, 0.0);
57 result.utilization.assign(ne, 0.0);
58 result.active.assign(ne, true);
59 return result; // trivially zero -> wrong
60 }
61
62 } // namespace truss2d
63
/opt/truss2d/include/truss2d/types.hpp
1 #ifndef TRUSS2D_TYPES_HPP
2 #define TRUSS2D_TYPES_HPP
3
4 #include <array>
5 #include <cstddef>
6
7 namespace truss2d {
8
9 /// A point / vector in the 2D plane.
10 struct Vec2 {
11 double x{0.0};
12 double y{0.0};
13 };
14
15 /// A structural node with a planar position. Each node owns two
16 /// translational degrees of freedom (DOFs): x then y.
17 struct Node {
18 Vec2 position{};
19 };
20
21 /// Linear-elastic axial material + section properties for a bar element.
22 /// `E` is Young's modulus [Pa], `A` is the cross-sectional area [m^2].
23 struct Section {
24 double E{0.0};
25 double A{0.0};
26 };
27
28 /// Force-transmission character of a member.
29 ///
30 /// Bar - a two-force member that resists BOTH tension and compression
31 /// (the classical bidirectional pin-jointed bar).
32 /// Cable - a slack-capable member that resists tension ONLY; it carries no
33 /// force and contributes no stiffness when it would otherwise be in
34 /// compression.
35 /// Strut - a contact-only member that resists compression ONLY; it carries
36 /// no force and contributes no stiffness when it would otherwise be
37 /// in tension (the ends separate).
38 enum class MemberKind { Bar, Cable, Strut };
39
40 /// A two-force member connecting node `n1` to node `n2`. Indices reference the
41 /// node array stored on the Model.
42 ///
43 /// `prestrain`, `alpha`, and `dT` describe installed strain sources. They
44 /// produce internal axial force even at zero nodal displacement; see the task
45 /// contract for the sign convention.
46 struct Element {
47 std::size_t n1{0};
48 std::size_t n2{0};
49 Section section{};
50 MemberKind kind{MemberKind::Bar};
51 double prestrain{0.0};
52 double alpha{0.0};
53 double dT{0.0};
54 };
55
56 /// A grounded linear spring support at a single node.
57 ///
58 /// The support resists displacement of node `node` ALONG the unit direction
59 /// `dir` with stiffness `kappa` [N/m]. It produces no resistance to motion
60 /// perpendicular to `dir` (a skew roller). Two springs on one node with
61 /// independent directions therefore behave like a 2D elastic support; a single
62 /// spring behaves like an inclined (skew) roller.
63 ///
64 /// The support may also SETTLE: its grounded end is held at a prescribed offset
65 /// `settlement` [m] measured along `dir`, so the spring's stored elongation is
66 /// the node's displacement along `dir` relative to that offset. With
67 /// `settlement = 0` this is an ordinary grounded spring.
68 ///
69 /// `dir` is stored normalized by the Model.
70 struct SpringSupport {
71 std::size_t node{0};
72 Vec2 dir{1.0, 0.0};
73 double kappa{0.0};
74 double settlement{0.0};
75 };
76
77 /// An applied nodal force. `dof` is the global DOF index
78 /// (2*node + component) and `value` is the force magnitude [N].
79 struct Load {
80 std::size_t dof{0};
81 double value{0.0};
82 };
83
84 /// The two global DOF indices owned by a node: {2*node, 2*node + 1}.
85 inline std::array<std::size_t, 2> node_dofs(std::size_t node) {
86 return {2 * node, 2 * node + 1};
87 }
88
89 } // namespace truss2d
90
91 #endif // TRUSS2D_TYPES_HPP
92
/opt/truss2d/include/truss2d/dense_matrix.hpp
1 #ifndef TRUSS2D_DENSE_MATRIX_HPP
2 #define TRUSS2D_DENSE_MATRIX_HPP
3
4 #include <cstddef>
5 #include <stdexcept>
6 #include <vector>
7
8 namespace truss2d {
9
10 /// A small, row-major dense matrix of doubles.
11 ///
12 /// This is part of the engine's linear-algebra support layer and is fully
13 /// implemented. The truss solver uses it to hold the global stiffness matrix
14 /// and to perform the linear solve via `solve_spd`.
15 class DenseMatrix {
16 public:
17 DenseMatrix() = default;
18
19 DenseMatrix(std::size_t rows, std::size_t cols)
20 : rows_(rows), cols_(cols), data_(rows * cols, 0.0) {}
21
22 std::size_t rows() const { return rows_; }
23 std::size_t cols() const { return cols_; }
24
25 double& operator()(std::size_t r, std::size_t c) {
26 return data_[r * cols_ + c];
27 }
28 double operator()(std::size_t r, std::size_t c) const {
29 return data_[r * cols_ + c];
30 }
31
32 void fill(double v) {
33 for (auto& x : data_) x = v;
34 }
35
36 private:
37 std::size_t rows_{0};
38 std::size_t cols_{0};
39 std::vector<double> data_;
40 };
41
42 /// Solve A x = b for a symmetric positive-definite matrix A using an
43 /// LDL^T (Cholesky-style) factorization with no pivoting.
44 ///
45 /// `A` is the n-by-n system matrix and `b` is the right-hand side of length n.
46 /// Returns the solution vector x of length n.
47 ///
48 /// Throws std::runtime_error if A is not square, if the dimensions are
49 /// inconsistent, or if A is detected to be singular / not positive-definite
50 /// (a zero or negative pivot is encountered). This factorization is
51 /// numerically suitable for the reduced stiffness matrix of a properly
52 /// constrained truss, which is SPD.
53 std::vector<double> solve_spd(const DenseMatrix& A, const std::vector<double>& b);
54
55 } // namespace truss2d
56
57 #endif // TRUSS2D_DENSE_MATRIX_HPP
58
/opt/truss2d/include/truss2d/model.hpp
1 #ifndef TRUSS2D_MODEL_HPP
2 #define TRUSS2D_MODEL_HPP
3
4 #include <cstddef>
5 #include <vector>
6
7 #include "truss2d/types.hpp"
8
9 namespace truss2d {
10
11 /// A 2D structural network of two-force members with spring (skew-roller)
12 /// supports and per-member tension/compression character and prestrain.
13 ///
14 /// The model is a plain data container plus light validation / derived-quantity
15 /// helpers. It is fully implemented; the numerical analysis lives in
16 /// StaticSolver (see solver.hpp).
17 class Model {
18 public:
19 std::size_t add_node(double x, double y);
20
21 /// Add a member between two existing nodes.
22 /// `kind` - Bar (default), Cable, or Strut (see types.hpp).
23 /// `prestrain` - mechanical installed axial strain (default 0).
24 /// `alpha,dT` - optional thermal expansion data (defaults 0).
25 /// Throws std::out_of_range if a node index is invalid.
26 std::size_t add_element(std::size_t n1, std::size_t n2, Section section,
27 MemberKind kind = MemberKind::Bar,
28 double prestrain = 0.0,
29 double alpha = 0.0,
30 double dT = 0.0);
31
32 /// Add a grounded spring support at `node` resisting displacement along
33 /// direction (dx, dy) with stiffness `kappa` [N/m]. The direction is
34 /// normalized internally. Throws std::out_of_range if `node` is invalid,
35 /// std::runtime_error if (dx, dy) is the zero vector or kappa <= 0.
36 /// `settlement` (optional, default 0) prescribes the spring's grounded-end
37 /// offset along (dx, dy); the spring force depends on the node displacement
38 /// along that direction RELATIVE to this offset.
39 void add_spring(std::size_t node, double dx, double dy, double kappa,
40 double settlement = 0.0);
41
42 /// Apply a force `value` [N] to a single global DOF (2*node + component).
43 /// Throws std::out_of_range if the DOF is invalid.
44 void add_load(std::size_t dof, double value);
45
46 const std::vector<Node>& nodes() const { return nodes_; }
47 const std::vector<Element>& elements() const { return elements_; }
48 const std::vector<SpringSupport>& springs() const { return springs_; }
49 const std::vector<Load>& loads() const { return loads_; }
50
51 std::size_t num_nodes() const { return nodes_.size(); }
52 std::size_t num_elements() const { return elements_.size(); }
53 std::size_t num_dofs() const { return 2 * nodes_.size(); }
54
55 /// Undeformed length of element `e` [m].
56 /// Throws std::out_of_range if `e` is invalid.
57 double element_length(std::size_t e) const;
58
59 /// Direction cosines (cos, sin) of element `e` measured from node n1
60 /// toward node n2, in the global frame.
61 /// Throws std::out_of_range if `e` is invalid, std::runtime_error if the
62 /// element has zero length.
63 Vec2 element_direction(std::size_t e) const;
64
65 /// Throws std::runtime_error if the model is structurally ill-formed
66 /// (no nodes, no elements, a zero-length element, a member with
67 /// non-positive E or A, or no spring supports at all).
68 void validate() const;
69
70 private:
71 std::vector<Node> nodes_;
72 std::vector<Element> elements_;
73 std::vector<SpringSupport> springs_;
74 std::vector<Load> loads_;
75 };
76
77 } // namespace truss2d
78
79 #endif // TRUSS2D_MODEL_HPP
80
/opt/truss2d/include/truss2d/types.hpp
1 #ifndef TRUSS2D_TYPES_HPP
2 #define TRUSS2D_TYPES_HPP
3
4 #include <array>
5 #include <cstddef>
6
7 namespace truss2d {
8
9 /// A point / vector in the 2D plane.
10 struct Vec2 {
11 double x{0.0};
12 double y{0.0};
13 };
14
15 /// A structural node with a planar position. Each node owns two
16 /// translational degrees of freedom (DOFs): x then y.
17 struct Node {
18 Vec2 position{};
19 };
20
21 /// Linear-elastic axial material + section properties for a bar element.
22 /// `E` is Young's modulus [Pa], `A` is the cross-sectional area [m^2].
23 struct Section {
24 double E{0.0};
25 double A{0.0};
26 };
27
28 /// Force-transmission character of a member.
29 ///
30 /// Bar - a two-force member that resists BOTH tension and compression
31 /// (the classical bidirectional pin-jointed bar).
32 /// Cable - a slack-capable member that resists tension ONLY; it carries no
33 /// force and contributes no stiffness when it would otherwise be in
34 /// compression.
35 /// Strut - a contact-only member that resists compression ONLY; it carries
36 /// no force and contributes no stiffness when it would otherwise be
37 /// in tension (the ends separate).
38 enum class MemberKind { Bar, Cable, Strut };
39
40 /// A two-force member connecting node `n1` to node `n2`. Indices reference the
41 /// node array stored on the Model.
42 ///
43 /// `prestrain`, `alpha`, and `dT` describe installed strain sources. They
44 /// produce internal axial force even at zero nodal displacement; see the task
45 /// contract for the sign convention.
46 struct Element {
47 std::size_t n1{0};
48 std::size_t n2{0};
49 Section section{};
50 MemberKind kind{MemberKind::Bar};
51 double prestrain{0.0};
52 double alpha{0.0};
53 double dT{0.0};
54 };
55
56 /// A grounded linear spring support at a single node.
57 ///
58 /// The support resists displacement of node `node` ALONG the unit direction
59 /// `dir` with stiffness `kappa` [N/m]. It produces no resistance to motion
60 /// perpendicular to `dir` (a skew roller). Two springs on one node with
61 /// independent directions therefore behave like a 2D elastic support; a single
62 /// spring behaves like an inclined (skew) roller.
63 ///
64 /// The support may also SETTLE: its grounded end is held at a prescribed offset
65 /// `settlement` [m] measured along `dir`, so the spring's stored elongation is
66 /// the node's displacement along `dir` relative to that offset. With
67 /// `settlement = 0` this is an ordinary grounded spring.
68 ///
69 /// `dir` is stored normalized by the Model.
70 struct SpringSupport {
71 std::size_t node{0};
72 Vec2 dir{1.0, 0.0};
73 double kappa{0.0};
74 double settlement{0.0};
75 };
76
77 /// An applied nodal force. `dof` is the global DOF index
78 /// (2*node + component) and `value` is the force magnitude [N].
79 struct Load {
80 std::size_t dof{0};
81 double value{0.0};
82 };
83
84 /// The two global DOF indices owned by a node: {2*node, 2*node + 1}.
85 inline std::array<std::size_t, 2> node_dofs(std::size_t node) {
86 return {2 * node, 2 * node + 1};
87 }
88
89 } // namespace truss2d
90
91 #endif // TRUSS2D_TYPES_HPP
92
/opt/truss2d/include/truss2d/dense_matrix.hpp
1 #ifndef TRUSS2D_DENSE_MATRIX_HPP
2 #define TRUSS2D_DENSE_MATRIX_HPP
3
4 #include <cstddef>
5 #include <stdexcept>
6 #include <vector>
7
8 namespace truss2d {
9
10 /// A small, row-major dense matrix of doubles.
11 ///
12 /// This is part of the engine's linear-algebra support layer and is fully
13 /// implemented. The truss solver uses it to hold the global stiffness matrix
14 /// and to perform the linear solve via `solve_spd`.
15 class DenseMatrix {
16 public:
17 DenseMatrix() = default;
18
19 DenseMatrix(std::size_t rows, std::size_t cols)
20 : rows_(rows), cols_(cols), data_(rows * cols, 0.0) {}
21
22 std::size_t rows() const { return rows_; }
23 std::size_t cols() const { return cols_; }
24
25 double& operator()(std::size_t r, std::size_t c) {
26 return data_[r * cols_ + c];
27 }
28 double operator()(std::size_t r, std::size_t c) const {
29 return data_[r * cols_ + c];
30 }
31
32 void fill(double v) {
33 for (auto& x : data_) x = v;
34 }
35
36 private:
37 std::size_t rows_{0};
38 std::size_t cols_{0};
39 std::vector<double> data_;
40 };
41
42 /// Solve A x = b for a symmetric positive-definite matrix A using an
43 /// LDL^T (Cholesky-style) factorization with no pivoting.
44 ///
45 /// `A` is the n-by-n system matrix and `b` is the right-hand side of length n.
46 /// Returns the solution vector x of length n.
47 ///
48 /// Throws std::runtime_error if A is not square, if the dimensions are
49 /// inconsistent, or if A is detected to be singular / not positive-definite
50 /// (a zero or negative pivot is encountered). This factorization is
51 /// numerically suitable for the reduced stiffness matrix of a properly
52 /// constrained truss, which is SPD.
53 std::vector<double> solve_spd(const DenseMatrix& A, const std::vector<double>& b);
54
55 } // namespace truss2d
56
57 #endif // TRUSS2D_DENSE_MATRIX_HPP
58
/opt/truss2d/src/model.cpp
1 #include "truss2d/model.hpp"
2
3 #include <cmath>
4 #include <stdexcept>
5
6 namespace truss2d {
7
8 std::size_t Model::add_node(double x, double y) {
9 nodes_.push_back(Node{Vec2{x, y}});
10 return nodes_.size() - 1;
11 }
12
13 std::size_t Model::add_element(std::size_t n1, std::size_t n2, Section section,
14 MemberKind kind, double prestrain,
15 double alpha, double dT) {
16 if (n1 >= nodes_.size() || n2 >= nodes_.size()) {
17 throw std::out_of_range("Model::add_element: node index out of range");
18 }
19 elements_.push_back(Element{n1, n2, section, kind, prestrain, alpha, dT});
20 return elements_.size() - 1;
21 }
22
23 void Model::add_spring(std::size_t node, double dx, double dy, double kappa,
24 double settlement) {
25 if (node >= nodes_.size()) {
26 throw std::out_of_range("Model::add_spring: node index out of range");
27 }
28 const double len = std::sqrt(dx * dx + dy * dy);
29 if (len <= 0.0) {
30 throw std::runtime_error("Model::add_spring: zero direction vector");
31 }
32 if (kappa <= 0.0) {
33 throw std::runtime_error("Model::add_spring: non-positive stiffness");
34 }
35 springs_.push_back(
36 SpringSupport{node, Vec2{dx / len, dy / len}, kappa, settlement});
37 }
38
39 void Model::add_load(std::size_t dof, double value) {
40 if (dof >= num_dofs()) {
41 throw std::out_of_range("Model::add_load: DOF out of range");
42 }
43 loads_.push_back(Load{dof, value});
44 }
45
46 double Model::element_length(std::size_t e) const {
47 if (e >= elements_.size()) {
48 throw std::out_of_range("Model::element_length: element index out of range");
49 }
50 const Element& el = elements_[e];
51 const Vec2& a = nodes_[el.n1].position;
52 const Vec2& b = nodes_[el.n2].position;
53 const double dx = b.x - a.x;
54 const double dy = b.y - a.y;
55 return std::sqrt(dx * dx + dy * dy);
56 }
57
58 Vec2 Model::element_direction(std::size_t e) const {
59 if (e >= elements_.size()) {
60 throw std::out_of_range("Model::element_direction: element index out of range");
61 }
62 const Element& el = elements_[e];
63 const Vec2& a = nodes_[el.n1].position;
64 const Vec2& b = nodes_[el.n2].position;
65 const double dx = b.x - a.x;
66 const double dy = b.y - a.y;
67 const double len = std::sqrt(dx * dx + dy * dy);
68 if (len <= 0.0) {
69 throw std::runtime_error("Model::element_direction: zero-length element");
70 }
71 return Vec2{dx / len, dy / len};
72 }
73
74 void Model::validate() const {
75 if (nodes_.empty()) {
76 throw std::runtime_error("Model::validate: model has no nodes");
77 }
78 if (elements_.empty()) {
79 throw std::runtime_error("Model::validate: model has no elements");
80 }
81 for (std::size_t e = 0; e < elements_.size(); ++e) {
82 const Element& el = elements_[e];
83 if (el.section.E <= 0.0 || el.section.A <= 0.0) {
84 throw std::runtime_error("Model::validate: element has non-positive E or A");
85 }
86 if (element_length(e) <= 0.0) {
87 throw std::runtime_error("Model::validate: element has zero length");
88 }
89 }
90 // Without at least one spring support the network is free to drift; full
91 // mechanism detection is left to the linear solve, which reports a singular
92 // system.
93 if (springs_.empty()) {
94 throw std::runtime_error(
95 "Model::validate: no spring supports (rigid-body motion not suppressed)");
96 }
97 }
98
99 } // namespace truss2d
100
/opt/truss2d/src/dense_matrix.cpp
1 #include "truss2d/dense_matrix.hpp"
2
3 #include <cmath>
4
5 namespace truss2d {
6
7 std::vector<double> solve_spd(const DenseMatrix& A, const std::vector<double>& b) {
8 const std::size_t n = A.rows();
9 if (A.cols() != n) {
10 throw std::runtime_error("solve_spd: matrix is not square");
11 }
12 if (b.size() != n) {
13 throw std::runtime_error("solve_spd: right-hand side size mismatch");
14 }
15
16 // LDL^T factorization (no pivoting): A = L D L^T with L unit-lower
17 // triangular and D diagonal. Suitable for the SPD reduced stiffness
18 // matrix of a well-constrained truss.
19 DenseMatrix L(n, n);
20 std::vector<double> D(n, 0.0);
21
22 for (std::size_t j = 0; j < n; ++j) {
23 double dj = A(j, j);
24 for (std::size_t k = 0; k < j; ++k) {
25 dj -= L(j, k) * L(j, k) * D[k];
26 }
27 // A genuine SPD matrix has strictly positive pivots. A tiny tolerance
28 // guards against round-off while still flagging singular systems.
29 if (dj <= 1e-12) {
30 throw std::runtime_error(
31 "solve_spd: matrix is singular or not positive-definite");
32 }
33 D[j] = dj;
34 L(j, j) = 1.0;
35 for (std::size_t i = j + 1; i < n; ++i) {
36 double s = A(i, j);
37 for (std::size_t k = 0; k < j; ++k) {
38 s -= L(i, k) * L(j, k) * D[k];
39 }
40 L(i, j) = s / dj;
41 }
42 }
43
44 // Forward solve L z = b.
45 std::vector<double> z(n, 0.0);
46 for (std::size_t i = 0; i < n; ++i) {
47 double s = b[i];
48 for (std::size_t k = 0; k < i; ++k) {
49 s -= L(i, k) * z[k];
50 }
51 z[i] = s;
52 }
53
54 // Diagonal solve D y = z.
55 std::vector<double> y(n, 0.0);
56 for (std::size_t i = 0; i < n; ++i) {
57 y[i] = z[i] / D[i];
58 }
59
60 // Back solve L^T x = y.
61 std::vector<double> x(n, 0.0);
62 for (std::size_t ii = 0; ii < n; ++ii) {
63 const std::size_t i = n - 1 - ii;
64 double s = y[i];
65 for (std::size_t k = i + 1; k < n; ++k) {
66 s -= L(k, i) * x[k];
67 }
68 x[i] = s;
69 }
70
71 return x;
72 }
73
74 } // namespace truss2d
75
/opt/truss2d/src/model.cpp
1 #include "truss2d/model.hpp"
2
3 #include <cmath>
4 #include <stdexcept>
5
6 namespace truss2d {
7
8 std::size_t Model::add_node(double x, double y) {
9 nodes_.push_back(Node{Vec2{x, y}});
10 return nodes_.size() - 1;
11 }
12
13 std::size_t Model::add_element(std::size_t n1, std::size_t n2, Section section,
14 MemberKind kind, double prestrain,
15 double alpha, double dT) {
16 if (n1 >= nodes_.size() || n2 >= nodes_.size()) {
17 throw std::out_of_range("Model::add_element: node index out of range");
18 }
19 elements_.push_back(Element{n1, n2, section, kind, prestrain, alpha, dT});
20 return elements_.size() - 1;
21 }
22
23 void Model::add_spring(std::size_t node, double dx, double dy, double kappa,
24 double settlement) {
25 if (node >= nodes_.size()) {
26 throw std::out_of_range("Model::add_spring: node index out of range");
27 }
28 const double len = std::sqrt(dx * dx + dy * dy);
29 if (len <= 0.0) {
30 throw std::runtime_error("Model::add_spring: zero direction vector");
31 }
32 if (kappa <= 0.0) {
33 throw std::runtime_error("Model::add_spring: non-positive stiffness");
34 }
35 springs_.push_back(
36 SpringSupport{node, Vec2{dx / len, dy / len}, kappa, settlement});
37 }
38
39 void Model::add_load(std::size_t dof, double value) {
40 if (dof >= num_dofs()) {
41 throw std::out_of_range("Model::add_load: DOF out of range");
42 }
43 loads_.push_back(Load{dof, value});
44 }
45
46 double Model::element_length(std::size_t e) const {
47 if (e >= elements_.size()) {
48 throw std::out_of_range("Model::element_length: element index out of range");
49 }
50 const Element& el = elements_[e];
51 const Vec2& a = nodes_[el.n1].position;
52 const Vec2& b = nodes_[el.n2].position;
53 const double dx = b.x - a.x;
54 const double dy = b.y - a.y;
55 return std::sqrt(dx * dx + dy * dy);
56 }
57
58 Vec2 Model::element_direction(std::size_t e) const {
59 if (e >= elements_.size()) {
60 throw std::out_of_range("Model::element_direction: element index out of range");
61 }
62 const Element& el = elements_[e];
63 const Vec2& a = nodes_[el.n1].position;
64 const Vec2& b = nodes_[el.n2].position;
65 const double dx = b.x - a.x;
66 const double dy = b.y - a.y;
67 const double len = std::sqrt(dx * dx + dy * dy);
68 if (len <= 0.0) {
69 throw std::runtime_error("Model::element_direction: zero-length element");
70 }
71 return Vec2{dx / len, dy / len};
72 }
73
74 void Model::validate() const {
75 if (nodes_.empty()) {
76 throw std::runtime_error("Model::validate: model has no nodes");
77 }
78 if (elements_.empty()) {
79 throw std::runtime_error("Model::validate: model has no elements");
80 }
81 for (std::size_t e = 0; e < elements_.size(); ++e) {
82 const Element& el = elements_[e];
83 if (el.section.E <= 0.0 || el.section.A <= 0.0) {
84 throw std::runtime_error("Model::validate: element has non-positive E or A");
85 }
86 if (element_length(e) <= 0.0) {
87 throw std::runtime_error("Model::validate: element has zero length");
88 }
89 }
90 // Without at least one spring support the network is free to drift; full
91 // mechanism detection is left to the linear solve, which reports a singular
92 // system.
93 if (springs_.empty()) {
94 throw std::runtime_error(
95 "Model::validate: no spring supports (rigid-body motion not suppressed)");
96 }
97 }
98
99 } // namespace truss2d
100
/opt/truss2d/src/dense_matrix.cpp
1 #include "truss2d/dense_matrix.hpp"
2
3 #include <cmath>
4
5 namespace truss2d {
6
7 std::vector<double> solve_spd(const DenseMatrix& A, const std::vector<double>& b) {
8 const std::size_t n = A.rows();
9 if (A.cols() != n) {
10 throw std::runtime_error("solve_spd: matrix is not square");
11 }
12 if (b.size() != n) {
13 throw std::runtime_error("solve_spd: right-hand side size mismatch");
14 }
15
16 // LDL^T factorization (no pivoting): A = L D L^T with L unit-lower
17 // triangular and D diagonal. Suitable for the SPD reduced stiffness
18 // matrix of a well-constrained truss.
19 DenseMatrix L(n, n);
20 std::vector<double> D(n, 0.0);
21
22 for (std::size_t j = 0; j < n; ++j) {
23 double dj = A(j, j);
24 for (std::size_t k = 0; k < j; ++k) {
25 dj -= L(j, k) * L(j, k) * D[k];
26 }
27 // A genuine SPD matrix has strictly positive pivots. A tiny tolerance
28 // guards against round-off while still flagging singular systems.
29 if (dj <= 1e-12) {
30 throw std::runtime_error(
31 "solve_spd: matrix is singular or not positive-definite");
32 }
33 D[j] = dj;
34 L(j, j) = 1.0;
35 for (std::size_t i = j + 1; i < n; ++i) {
36 double s = A(i, j);
37 for (std::size_t k = 0; k < j; ++k) {
38 s -= L(i, k) * L(j, k) * D[k];
39 }
40 L(i, j) = s / dj;
41 }
42 }
43
44 // Forward solve L z = b.
45 std::vector<double> z(n, 0.0);
46 for (std::size_t i = 0; i < n; ++i) {
47 double s = b[i];
48 for (std::size_t k = 0; k < i; ++k) {
49 s -= L(i, k) * z[k];
50 }
51 z[i] = s;
52 }
53
54 // Diagonal solve D y = z.
55 std::vector<double> y(n, 0.0);
56 for (std::size_t i = 0; i < n; ++i) {
57 y[i] = z[i] / D[i];
58 }
59
60 // Back solve L^T x = y.
61 std::vector<double> x(n, 0.0);
62 for (std::size_t ii = 0; ii < n; ++ii) {
63 const std::size_t i = n - 1 - ii;
64 double s = y[i];
65 for (std::size_t k = i + 1; k < n; ++k) {
66 s -= L(k, i) * x[k];
67 }
68 x[i] = s;
69 }
70
71 return x;
72 }
73
74 } // namespace truss2d
75
/opt/truss2d/instruction.md
File does not exist. Note: your current working directory is /opt/truss2d. [error] tool reported failure
List project tree
ls -R /opt/truss2d --ignore=build
/opt/truss2d: CMakeLists.txt examples include src tests /opt/truss2d/examples: two_bar.truss /opt/truss2d/include: truss2d /opt/truss2d/include/truss2d: dense_matrix.hpp io.hpp model.hpp solver.hpp types.hpp /opt/truss2d/src: dense_matrix.cpp io.cpp main.cpp model.cpp solver.cpp /opt/truss2d/tests: CMakeLists.txt test_smoke.cpp test_util.hpp [stdout] /opt/truss2d: CMakeLists.txt examples include src tests /opt/truss2d/examples: two_bar.truss /opt/truss2d/include: truss2d /opt/truss2d/include/truss2d: dense_matrix.hpp io.hpp model.hpp solver.hpp types.hpp /opt/truss2d/src: dense_matrix.cpp io.cpp main.cpp model.cpp solver.cpp /opt/truss2d/tests: CMakeLists.txt test_smoke.cpp test_util.hpp
ls -la /opt/truss2d && echo "---" && find /opt/truss2d -name "*.md" && echo "---tests---" && ls -R /opt/truss2d/tests 2>/dev/null && echo "---apps---" && ls -R /opt/truss2d/apps /opt/truss2d/src /opt/truss2d/examples 2>/dev/null
Exit code 2 total 4 drwxr-xr-x 6 root root 128 Jun 21 09:20 . drwxr-xr-x 3 root root 60 Jun 20 18:52 .. -rw-r--r-- 1 root root 807 Jun 20 18:51 CMakeLists.txt drwxr-xr-x 2 root root 60 Jun 20 18:51 examples drwxr-xr-x 3 root root 60 Jun 20 18:51 include drwxr-xr-x 2 root root 128 Jun 21 09:13 src drwxr-xr-x 2 root root 128 Jun 20 18:51 tests --- ---tests--- /opt/truss2d/tests: CMakeLists.txt test_smoke.cpp test_util.hpp ---apps--- /opt/truss2d/examples: two_bar.truss /opt/truss2d/src: dense_matrix.cpp io.cpp main.cpp model.cpp solver.cpp [error] tool reported failure
/opt/truss2d/tests/test_smoke.cpp
1 // Trivial smoke test for the truss2d engine.
2 //
3 // This checks only that the engine links and that the data model + result
4 // shapes are self-consistent. It deliberately makes NO numerical-correctness
5 // assertions, so it passes against the shipped stub solver. The real
6 // correctness / invariant / discriminator suite is provided separately by the
7 // grader.
8
9 #include <cstddef>
10 #include <sstream>
11
12 #include "truss2d/io.hpp"
13 #include "truss2d/model.hpp"
14 #include "truss2d/solver.hpp"
15 #include "test_util.hpp"
16
17 using truss2d::Model;
18 using truss2d::Section;
19 using truss2d::MemberKind;
20 using truss2d::SolveResult;
21 using truss2d::StaticSolver;
22
23 TRUSS2D_TEST("model_construction_shapes") {
24 Model m;
25 m.add_node(0.0, 0.0);
26 m.add_node(2.0, 0.0);
27 m.add_element(0, 1, Section{2.0e11, 1.0e-3}, MemberKind::Bar, 0.0);
28 m.add_spring(0, 1.0, 0.0, 1.0e13);
29 m.add_spring(0, 0.0, 1.0, 1.0e13);
30 m.add_spring(1, 0.0, 1.0, 1.0e13);
31 m.add_load(2, 1.0e3);
32
33 truss2d_test::expect(m.num_nodes() == 2, "two nodes");
34 truss2d_test::expect(m.num_dofs() == 4, "four DOFs");
35 truss2d_test::expect(m.num_elements() == 1, "one element");
36 truss2d_test::expect(m.springs().size() == 3, "three springs");
37 }
38
39 TRUSS2D_TEST("result_vectors_have_expected_lengths") {
40 Model m;
41 m.add_node(0.0, 0.0);
42 m.add_node(3.0, 0.0);
43 m.add_node(3.0, 4.0);
44 m.add_element(0, 1, Section{2.0e11, 1.0e-3});
45 m.add_element(1, 2, Section{2.0e11, 1.0e-3});
46 m.add_element(0, 2, Section{2.0e11, 1.0e-3});
47 m.add_spring(0, 1.0, 0.0, 1.0e13);
48 m.add_spring(0, 0.0, 1.0, 1.0e13);
49 m.add_spring(1, 0.0, 1.0, 1.0e13);
50
51 SolveResult r = StaticSolver(m).solve();
52 truss2d_test::expect(r.displacements.size() == m.num_dofs(),
53 "displacements length == num_dofs");
54 truss2d_test::expect(r.reactions.size() == m.num_dofs(),
55 "reactions length == num_dofs");
56 truss2d_test::expect(r.axial_forces.size() == m.num_elements(),
57 "axial_forces length == num_elements");
58 truss2d_test::expect(r.utilization.size() == m.num_elements(),
59 "utilization length == num_elements");
60 truss2d_test::expect(r.active.size() == m.num_elements(),
61 "active length == num_elements");
62 }
63
64 TRUSS2D_TEST("parser_reads_records") {
65 const char* text =
66 "node 0 0\n"
67 "node 2 0\n"
68 "elem 0 1 2.0e11 1.0e-3 1 -5.0e-4 1.2e-5 -25\n"
69 "spring 0 1 0 1.0e13\n"
70 "spring 0 0 1 1.0e13\n"
71 "spring 1 0 1 1.0e13\n"
72 "load 1 1 -1.0e3\n";
73 std::istringstream in(text);
74 Model m = truss2d::parse_model(in);
75 truss2d_test::expect(m.num_nodes() == 2, "parsed two nodes");
76 truss2d_test::expect(m.num_elements() == 1, "parsed one element");
77 truss2d_test::expect(m.elements()[0].kind == MemberKind::Cable,
78 "parsed cable kind");
79 truss2d_test::expect_near(m.elements()[0].alpha, 1.2e-5, 1e-15,
80 "parsed alpha");
81 truss2d_test::expect_near(m.elements()[0].dT, -25.0, 1e-15,
82 "parsed dT");
83 truss2d_test::expect(m.springs().size() == 3, "parsed three springs");
84 }
85
86 TRUSS2D_TEST_MAIN()
87
/opt/truss2d/tests/test_util.hpp
1 #ifndef TRUSS2D_TEST_UTIL_HPP
2 #define TRUSS2D_TEST_UTIL_HPP
3
4 // Minimal dependency-free test harness for the truss2d suite.
5 //
6 // Each test executable defines test functions and registers them with
7 // TRUSS2D_TEST(name) { ... }. main() (provided by test_main.hpp via
8 // TRUSS2D_TEST_MAIN) runs every registered test and returns non-zero if any
9 // assertion fails, which CTest interprets as a failure.
10
11 #include <cmath>
12 #include <cstddef>
13 #include <functional>
14 #include <iostream>
15 #include <string>
16 #include <vector>
17
18 namespace truss2d_test {
19
20 struct Case {
21 std::string name;
22 std::function<void()> fn;
23 };
24
25 inline std::vector<Case>& registry() {
26 static std::vector<Case> cases;
27 return cases;
28 }
29
30 struct Registrar {
31 Registrar(const std::string& name, std::function<void()> fn) {
32 registry().push_back(Case{name, std::move(fn)});
33 }
34 };
35
36 // Thrown by an assertion failure to abort the current test case.
37 struct AssertionError {
38 std::string message;
39 };
40
41 inline void expect(bool cond, const std::string& what) {
42 if (!cond) {
43 throw AssertionError{what};
44 }
45 }
46
47 inline void expect_near(double got, double want, double tol, const std::string& what) {
48 if (std::fabs(got - want) > tol) {
49 throw AssertionError{what + " (got " + std::to_string(got) +
50 ", want " + std::to_string(want) +
51 ", tol " + std::to_string(tol) + ")"};
52 }
53 }
54
55 inline int run_all() {
56 int failures = 0;
57 for (const auto& c : registry()) {
58 try {
59 c.fn();
60 std::cout << "[ PASS ] " << c.name << "\n";
61 } catch (const AssertionError& e) {
62 std::cout << "[ FAIL ] " << c.name << ": " << e.message << "\n";
63 ++failures;
64 } catch (const std::exception& e) {
65 std::cout << "[ FAIL ] " << c.name
66 << ": unexpected exception: " << e.what() << "\n";
67 ++failures;
68 } catch (...) {
69 std::cout << "[ FAIL ] " << c.name << ": unknown exception\n";
70 ++failures;
71 }
72 }
73 std::cout << "----\n"
74 << (registry().size() - failures) << "/" << registry().size()
75 << " tests passed\n";
76 return failures == 0 ? 0 : 1;
77 }
78
79 } // namespace truss2d_test
80
81 #define TRUSS2D_CONCAT_INNER(a, b) a##b
82 #define TRUSS2D_CONCAT(a, b) TRUSS2D_CONCAT_INNER(a, b)
83
84 #define TRUSS2D_TEST(NAME) \
85 static void TRUSS2D_CONCAT(truss2d_test_fn_, __LINE__)(); \
86 static ::truss2d_test::Registrar TRUSS2D_CONCAT(truss2d_test_reg_, \
87 __LINE__){ \
88 NAME, &TRUSS2D_CONCAT(truss2d_test_fn_, __LINE__)}; \
89 static void TRUSS2D_CONCAT(truss2d_test_fn_, __LINE__)()
90
91 #define TRUSS2D_TEST_MAIN() \
92 int main() { return ::truss2d_test::run_all(); }
93
94 #endif // TRUSS2D_TEST_UTIL_HPP
95
/opt/truss2d/src/main.cpp
1 #include <cstdlib>
2 #include <iostream>
3
4 #include "truss2d/io.hpp"
5 #include "truss2d/model.hpp"
6 #include "truss2d/solver.hpp"
7
8 // Command-line demo driver: read a .truss model, solve it, print a report.
9 // truss2d_cli <model.truss>
10 int main(int argc, char** argv) {
11 if (argc != 2) {
12 std::cerr << "usage: " << argv[0] << " <model.truss>\n";
13 return 2;
14 }
15 try {
16 truss2d::Model model = truss2d::parse_model_file(argv[1]);
17 truss2d::StaticSolver solver(model);
18 truss2d::SolveResult result = solver.solve();
19 truss2d::write_report(std::cout, model, result);
20 } catch (const std::exception& ex) {
21 std::cerr << "error: " << ex.what() << "\n";
22 return 1;
23 }
24 return 0;
25 }
26
/opt/truss2d/tests/test_smoke.cpp
1 // Trivial smoke test for the truss2d engine.
2 //
3 // This checks only that the engine links and that the data model + result
4 // shapes are self-consistent. It deliberately makes NO numerical-correctness
5 // assertions, so it passes against the shipped stub solver. The real
6 // correctness / invariant / discriminator suite is provided separately by the
7 // grader.
8
9 #include <cstddef>
10 #include <sstream>
11
12 #include "truss2d/io.hpp"
13 #include "truss2d/model.hpp"
14 #include "truss2d/solver.hpp"
15 #include "test_util.hpp"
16
17 using truss2d::Model;
18 using truss2d::Section;
19 using truss2d::MemberKind;
20 using truss2d::SolveResult;
21 using truss2d::StaticSolver;
22
23 TRUSS2D_TEST("model_construction_shapes") {
24 Model m;
25 m.add_node(0.0, 0.0);
26 m.add_node(2.0, 0.0);
27 m.add_element(0, 1, Section{2.0e11, 1.0e-3}, MemberKind::Bar, 0.0);
28 m.add_spring(0, 1.0, 0.0, 1.0e13);
29 m.add_spring(0, 0.0, 1.0, 1.0e13);
30 m.add_spring(1, 0.0, 1.0, 1.0e13);
31 m.add_load(2, 1.0e3);
32
33 truss2d_test::expect(m.num_nodes() == 2, "two nodes");
34 truss2d_test::expect(m.num_dofs() == 4, "four DOFs");
35 truss2d_test::expect(m.num_elements() == 1, "one element");
36 truss2d_test::expect(m.springs().size() == 3, "three springs");
37 }
38
39 TRUSS2D_TEST("result_vectors_have_expected_lengths") {
40 Model m;
41 m.add_node(0.0, 0.0);
42 m.add_node(3.0, 0.0);
43 m.add_node(3.0, 4.0);
44 m.add_element(0, 1, Section{2.0e11, 1.0e-3});
45 m.add_element(1, 2, Section{2.0e11, 1.0e-3});
46 m.add_element(0, 2, Section{2.0e11, 1.0e-3});
47 m.add_spring(0, 1.0, 0.0, 1.0e13);
48 m.add_spring(0, 0.0, 1.0, 1.0e13);
49 m.add_spring(1, 0.0, 1.0, 1.0e13);
50
51 SolveResult r = StaticSolver(m).solve();
52 truss2d_test::expect(r.displacements.size() == m.num_dofs(),
53 "displacements length == num_dofs");
54 truss2d_test::expect(r.reactions.size() == m.num_dofs(),
55 "reactions length == num_dofs");
56 truss2d_test::expect(r.axial_forces.size() == m.num_elements(),
57 "axial_forces length == num_elements");
58 truss2d_test::expect(r.utilization.size() == m.num_elements(),
59 "utilization length == num_elements");
60 truss2d_test::expect(r.active.size() == m.num_elements(),
61 "active length == num_elements");
62 }
63
64 TRUSS2D_TEST("parser_reads_records") {
65 const char* text =
66 "node 0 0\n"
67 "node 2 0\n"
68 "elem 0 1 2.0e11 1.0e-3 1 -5.0e-4 1.2e-5 -25\n"
69 "spring 0 1 0 1.0e13\n"
70 "spring 0 0 1 1.0e13\n"
71 "spring 1 0 1 1.0e13\n"
72 "load 1 1 -1.0e3\n";
73 std::istringstream in(text);
74 Model m = truss2d::parse_model(in);
75 truss2d_test::expect(m.num_nodes() == 2, "parsed two nodes");
76 truss2d_test::expect(m.num_elements() == 1, "parsed one element");
77 truss2d_test::expect(m.elements()[0].kind == MemberKind::Cable,
78 "parsed cable kind");
79 truss2d_test::expect_near(m.elements()[0].alpha, 1.2e-5, 1e-15,
80 "parsed alpha");
81 truss2d_test::expect_near(m.elements()[0].dT, -25.0, 1e-15,
82 "parsed dT");
83 truss2d_test::expect(m.springs().size() == 3, "parsed three springs");
84 }
85
86 TRUSS2D_TEST_MAIN()
87
/opt/truss2d/src/io.cpp
1 #include "truss2d/io.hpp"
2
3 #include <fstream>
4 #include <iomanip>
5 #include <iostream>
6 #include <sstream>
7 #include <stdexcept>
8 #include <string>
9
10 namespace truss2d {
11
12 namespace {
13
14 [[noreturn]] void fail(std::size_t line_no, const std::string& msg) {
15 std::ostringstream oss;
16 oss << "parse_model: line " << line_no << ": " << msg;
17 throw std::runtime_error(oss.str());
18 }
19
20 const char* kind_name(MemberKind k) {
21 switch (k) {
22 case MemberKind::Bar: return "bar";
23 case MemberKind::Cable: return "cable";
24 case MemberKind::Strut: return "strut";
25 }
26 return "?";
27 }
28
29 } // namespace
30
31 Model parse_model(std::istream& in) {
32 Model model;
33 std::string line;
34 std::size_t line_no = 0;
35
36 while (std::getline(in, line)) {
37 ++line_no;
38 const auto hash = line.find('#');
39 if (hash != std::string::npos) {
40 line.erase(hash);
41 }
42 std::istringstream ls(line);
43 std::string tag;
44 if (!(ls >> tag)) {
45 continue; // blank / comment-only line
46 }
47
48 if (tag == "node") {
49 double x, y;
50 if (!(ls >> x >> y)) fail(line_no, "expected: node <x> <y>");
51 model.add_node(x, y);
52 } else if (tag == "elem") {
53 std::size_t n1, n2;
54 double E, A;
55 if (!(ls >> n1 >> n2 >> E >> A)) {
56 fail(line_no,
57 "expected: elem <n1> <n2> <E> <A> [kind] [prestrain] [alpha] [dT]");
58 }
59 int kind_code = 0;
60 double prestrain = 0.0;
61 double alpha = 0.0;
62 double dT = 0.0;
63 ls >> kind_code; // optional
64 ls >> prestrain; // optional
65 ls >> alpha; // optional
66 ls >> dT; // optional
67 MemberKind kind = MemberKind::Bar;
68 if (kind_code == 1) kind = MemberKind::Cable;
69 else if (kind_code == 2) kind = MemberKind::Strut;
70 else if (kind_code != 0) fail(line_no, "kind must be 0/1/2 (bar/cable/strut)");
71 try {
72 model.add_element(n1, n2, Section{E, A}, kind, prestrain, alpha, dT);
73 } catch (const std::exception& ex) {
74 fail(line_no, ex.what());
75 }
76 } else if (tag == "spring") {
77 std::size_t node;
78 double dx, dy, kappa;
79 if (!(ls >> node >> dx >> dy >> kappa)) {
80 fail(line_no, "expected: spring <node> <dx> <dy> <kappa>");
81 }
82 try {
83 model.add_spring(node, dx, dy, kappa);
84 } catch (const std::exception& ex) {
85 fail(line_no, ex.what());
86 }
87 } else if (tag == "load") {
88 std::size_t node, comp;
89 double value;
90 if (!(ls >> node >> comp >> value)) {
91 fail(line_no, "expected: load <node> <comp> <value>");
92 }
93 if (comp > 1) fail(line_no, "component must be 0 (x) or 1 (y)");
94 try {
95 model.add_load(2 * node + comp, value);
96 } catch (const std::exception& ex) {
97 fail(line_no, ex.what());
98 }
99 } else {
100 fail(line_no, "unknown record '" + tag + "'");
101 }
102 }
103
104 return model;
105 }
106
107 Model parse_model_file(const std::string& path) {
108 std::ifstream in(path);
109 if (!in) {
110 throw std::runtime_error("parse_model_file: cannot open '" + path + "'");
111 }
112 return parse_model(in);
113 }
114
115 void write_report(std::ostream& out, const Model& model, const SolveResult& result) {
116 out << std::scientific << std::setprecision(6);
117 out << "# Truss2D static analysis report\n";
118 out << "nodes=" << model.num_nodes()
119 << " elements=" << model.num_elements()
120 << " dofs=" << model.num_dofs() << "\n\n";
121
122 out << "[displacements] (node, ux, uy)\n";
123 for (std::size_t i = 0; i < model.num_nodes(); ++i) {
124 out << i << " " << result.displacements[2 * i] << " "
125 << result.displacements[2 * i + 1] << "\n";
126 }
127
128 out << "\n[reactions] (node, Rx, Ry)\n";
129 for (std::size_t i = 0; i < model.num_nodes(); ++i) {
130 out << i << " " << result.reactions[2 * i] << " "
131 << result.reactions[2 * i + 1] << "\n";
132 }
133
134 out << "\n[axial_forces] (element, kind, active, N [+tension], utilization)\n";
135 for (std::size_t e = 0; e < model.num_elements(); ++e) {
136 const bool act = e < result.active.size() ? result.active[e] : true;
137 const double util = e < result.utilization.size() ? result.utilization[e] : 0.0;
138 out << e << " " << kind_name(model.elements()[e].kind) << " "
139 << (act ? 1 : 0) << " " << result.axial_forces[e] << " "
140 << util << "\n";
141 }
142 }
143
144 } // namespace truss2d
145
/opt/truss2d/src/main.cpp
1 #include <cstdlib>
2 #include <iostream>
3
4 #include "truss2d/io.hpp"
5 #include "truss2d/model.hpp"
6 #include "truss2d/solver.hpp"
7
8 // Command-line demo driver: read a .truss model, solve it, print a report.
9 // truss2d_cli <model.truss>
10 int main(int argc, char** argv) {
11 if (argc != 2) {
12 std::cerr << "usage: " << argv[0] << " <model.truss>\n";
13 return 2;
14 }
15 try {
16 truss2d::Model model = truss2d::parse_model_file(argv[1]);
17 truss2d::StaticSolver solver(model);
18 truss2d::SolveResult result = solver.solve();
19 truss2d::write_report(std::cout, model, result);
20 } catch (const std::exception& ex) {
21 std::cerr << "error: " << ex.what() << "\n";
22 return 1;
23 }
24 return 0;
25 }
26
/opt/truss2d/include/truss2d/io.hpp
1 #ifndef TRUSS2D_IO_HPP
2 #define TRUSS2D_IO_HPP
3
4 #include <iosfwd>
5 #include <string>
6
7 #include "truss2d/model.hpp"
8 #include "truss2d/solver.hpp"
9
10 namespace truss2d {
11
12 /// Parse a model from the engine's plain-text `.truss` format.
13 ///
14 /// The format is line-oriented. Blank lines and lines beginning with '#'
15 /// are ignored. Recognized records:
16 /// node <x> <y>
17 /// elem <n1> <n2> <E> <A> [<kind>] [<prestrain>] [<alpha>] [<dT>]
18 /// kind: 0=bar (default), 1=cable, 2=strut; trailing values default 0
19 /// spring <node> <dx> <dy> <kappa> # grounded spring along (dx,dy)
20 /// load <node> <comp> <value> # comp: 0 = x, 1 = y
21 ///
22 /// Node and element indices are assigned in declaration order starting at 0.
23 /// Throws std::runtime_error on malformed input.
24 Model parse_model(std::istream& in);
25
26 /// Convenience overload that parses from a file path.
27 Model parse_model_file(const std::string& path);
28
29 /// Write a human-readable analysis report (displacements, reactions, axial
30 /// forces, utilization, active flags) to `out`. Fully implemented; used by the
31 /// demo executable.
32 void write_report(std::ostream& out, const Model& model, const SolveResult& result);
33
34 } // namespace truss2d
35
36 #endif // TRUSS2D_IO_HPP
37
/opt/truss2d/CMakeLists.txt
1 cmake_minimum_required(VERSION 3.16)
2 project(truss2d LANGUAGES CXX)
3
4 set(CMAKE_CXX_STANDARD 17)
5 set(CMAKE_CXX_STANDARD_REQUIRED ON)
6 set(CMAKE_CXX_EXTENSIONS OFF)
7
8 if(NOT CMAKE_BUILD_TYPE)
9 set(CMAKE_BUILD_TYPE Release)
10 endif()
11
12 add_compile_options(-Wall -Wextra)
13
14 # --- Core library -----------------------------------------------------------
15 add_library(truss2d
16 src/dense_matrix.cpp
17 src/model.cpp
18 src/solver.cpp
19 src/io.cpp
20 )
21 target_include_directories(truss2d PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}/include)
22
23 # --- Demo CLI ---------------------------------------------------------------
24 add_executable(truss2d_cli src/main.cpp)
25 target_link_libraries(truss2d_cli PRIVATE truss2d)
26
27 # --- Tests ------------------------------------------------------------------
28 enable_testing()
29 add_subdirectory(tests)
30
/opt/truss2d/src/io.cpp
1 #include "truss2d/io.hpp"
2
3 #include <fstream>
4 #include <iomanip>
5 #include <iostream>
6 #include <sstream>
7 #include <stdexcept>
8 #include <string>
9
10 namespace truss2d {
11
12 namespace {
13
14 [[noreturn]] void fail(std::size_t line_no, const std::string& msg) {
15 std::ostringstream oss;
16 oss << "parse_model: line " << line_no << ": " << msg;
17 throw std::runtime_error(oss.str());
18 }
19
20 const char* kind_name(MemberKind k) {
21 switch (k) {
22 case MemberKind::Bar: return "bar";
23 case MemberKind::Cable: return "cable";
24 case MemberKind::Strut: return "strut";
25 }
26 return "?";
27 }
28
29 } // namespace
30
31 Model parse_model(std::istream& in) {
32 Model model;
33 std::string line;
34 std::size_t line_no = 0;
35
36 while (std::getline(in, line)) {
37 ++line_no;
38 const auto hash = line.find('#');
39 if (hash != std::string::npos) {
40 line.erase(hash);
41 }
42 std::istringstream ls(line);
43 std::string tag;
44 if (!(ls >> tag)) {
45 continue; // blank / comment-only line
46 }
47
48 if (tag == "node") {
49 double x, y;
50 if (!(ls >> x >> y)) fail(line_no, "expected: node <x> <y>");
51 model.add_node(x, y);
52 } else if (tag == "elem") {
53 std::size_t n1, n2;
54 double E, A;
55 if (!(ls >> n1 >> n2 >> E >> A)) {
56 fail(line_no,
57 "expected: elem <n1> <n2> <E> <A> [kind] [prestrain] [alpha] [dT]");
58 }
59 int kind_code = 0;
60 double prestrain = 0.0;
61 double alpha = 0.0;
62 double dT = 0.0;
63 ls >> kind_code; // optional
64 ls >> prestrain; // optional
65 ls >> alpha; // optional
66 ls >> dT; // optional
67 MemberKind kind = MemberKind::Bar;
68 if (kind_code == 1) kind = MemberKind::Cable;
69 else if (kind_code == 2) kind = MemberKind::Strut;
70 else if (kind_code != 0) fail(line_no, "kind must be 0/1/2 (bar/cable/strut)");
71 try {
72 model.add_element(n1, n2, Section{E, A}, kind, prestrain, alpha, dT);
73 } catch (const std::exception& ex) {
74 fail(line_no, ex.what());
75 }
76 } else if (tag == "spring") {
77 std::size_t node;
78 double dx, dy, kappa;
79 if (!(ls >> node >> dx >> dy >> kappa)) {
80 fail(line_no, "expected: spring <node> <dx> <dy> <kappa>");
81 }
82 try {
83 model.add_spring(node, dx, dy, kappa);
84 } catch (const std::exception& ex) {
85 fail(line_no, ex.what());
86 }
87 } else if (tag == "load") {
88 std::size_t node, comp;
89 double value;
90 if (!(ls >> node >> comp >> value)) {
91 fail(line_no, "expected: load <node> <comp> <value>");
92 }
93 if (comp > 1) fail(line_no, "component must be 0 (x) or 1 (y)");
94 try {
95 model.add_load(2 * node + comp, value);
96 } catch (const std::exception& ex) {
97 fail(line_no, ex.what());
98 }
99 } else {
100 fail(line_no, "unknown record '" + tag + "'");
101 }
102 }
103
104 return model;
105 }
106
107 Model parse_model_file(const std::string& path) {
108 std::ifstream in(path);
109 if (!in) {
110 throw std::runtime_error("parse_model_file: cannot open '" + path + "'");
111 }
112 return parse_model(in);
113 }
114
115 void write_report(std::ostream& out, const Model& model, const SolveResult& result) {
116 out << std::scientific << std::setprecision(6);
117 out << "# Truss2D static analysis report\n";
118 out << "nodes=" << model.num_nodes()
119 << " elements=" << model.num_elements()
120 << " dofs=" << model.num_dofs() << "\n\n";
121
122 out << "[displacements] (node, ux, uy)\n";
123 for (std::size_t i = 0; i < model.num_nodes(); ++i) {
124 out << i << " " << result.displacements[2 * i] << " "
125 << result.displacements[2 * i + 1] << "\n";
126 }
127
128 out << "\n[reactions] (node, Rx, Ry)\n";
129 for (std::size_t i = 0; i < model.num_nodes(); ++i) {
130 out << i << " " << result.reactions[2 * i] << " "
131 << result.reactions[2 * i + 1] << "\n";
132 }
133
134 out << "\n[axial_forces] (element, kind, active, N [+tension], utilization)\n";
135 for (std::size_t e = 0; e < model.num_elements(); ++e) {
136 const bool act = e < result.active.size() ? result.active[e] : true;
137 const double util = e < result.utilization.size() ? result.utilization[e] : 0.0;
138 out << e << " " << kind_name(model.elements()[e].kind) << " "
139 << (act ? 1 : 0) << " " << result.axial_forces[e] << " "
140 << util << "\n";
141 }
142 }
143
144 } // namespace truss2d
145
#include "truss2d/solver.hpp"
#include <array>
#include <cmath>
#include <cstddef>
#include <stdexcept>
#include <vector>
#include "truss2d/dense_matrix.hpp"
#include "truss2d/types.hpp"
namespace truss2d {
// ============================================================================
// Static solver for a 2D network of two-force members (bars/cables/struts) on
// finite spring supports, with installed (mechanical + thermal) strain and the
// small-displacement initial-stress (geometric) stiffening of taut members.
//
// Kinematics (small displacement). For element e with direction cosines
// (c, s) from n1 -> n2 and undeformed length L0, the axial extension is the
// first-order change of end separation projected onto the member axis:
//
// delta = c*(u2x - u1x) + s*(u2y - u1y) = b . u_local,
// b = (-c, -s, c, s).
//
// Axial force (tension positive). The installed free strain is
// eps0 = prestrain + alpha*dT
// (a member that is too long / heated has eps0 > 0 and is compressive when
// restrained). The elastic strain is the actual extension minus the free
// strain, so
// N = E*A * (delta / L0 - eps0).
//
// Governing system. For a given active set,
// ( K_elastic(active) + K_spring + K_geom(N) ) u = f_load + f_settle
// + f_prestrain(active),
// where the geometric stiffness of an active member is N/L0 times the outer
// product of the transverse kinematic vector b_t = (s, -c, -s, c), so it acts
// only on relative transverse motion, vanishes for rigid translation and for
// motion along the axis, stiffens under tension and softens under compression.
// K_geom depends on N which depends on u, so solve() iterates the member forces
// to a self-consistent state. The matrix helpers below report the ordinary
// elastic behaviour only; the geometric term lives entirely inside solve().
// ============================================================================
namespace {
constexpr double kAllowableStress = 2.5e8; // [Pa]
// Global DOF indices for an element's local DOFs {n1.x, n1.y, n2.x, n2.y}.
std::array<std::size_t, 4> element_dofs(const Element& el) {
return {2 * el.n1, 2 * el.n1 + 1, 2 * el.n2, 2 * el.n2 + 1};
}
} // namespace
DenseMatrix StaticSolver::element_stiffness_global(std::size_t e) const {
const Element& el = model_.elements()[e];
const double L = model_.element_length(e);
const Vec2 dir = model_.element_direction(e);
const double c = dir.x;
const double s = dir.y;
const double k = el.section.E * el.section.A / L;
// K_e = (E A / L) * b^T b, with b = (-c, -s, c, s).
const double b[4] = {-c, -s, c, s};
DenseMatrix Ke(4, 4);
for (std::size_t i = 0; i < 4; ++i) {
for (std::size_t j = 0; j < 4; ++j) {
Ke(i, j) = k * b[i] * b[j];
}
}
return Ke;
}
DenseMatrix StaticSolver::assemble_system(const std::vector<bool>& active) const {
const std::size_t n = model_.num_dofs();
const auto& elems = model_.elements();
DenseMatrix K(n, n);
// Elastic stiffness of every participating member.
for (std::size_t e = 0; e < elems.size(); ++e) {
if (!active[e]) continue;
const DenseMatrix Ke = element_stiffness_global(e);
const auto g = element_dofs(elems[e]);
for (std::size_t i = 0; i < 4; ++i) {
for (std::size_t j = 0; j < 4; ++j) {
K(g[i], g[j]) += Ke(i, j);
}
}
}
// Spring supports: kappa * dir dir^T in the supported node's 2x2 block.
for (const SpringSupport& sp : model_.springs()) {
const std::size_t gx = 2 * sp.node;
const std::size_t gy = 2 * sp.node + 1;
const double kx = sp.kappa * sp.dir.x;
const double ky = sp.kappa * sp.dir.y;
K(gx, gx) += kx * sp.dir.x;
K(gx, gy) += kx * sp.dir.y;
K(gy, gx) += ky * sp.dir.x;
K(gy, gy) += ky * sp.dir.y;
}
return K;
}
std::vector<double>
StaticSolver::prestrain_force(const std::vector<bool>& active) const {
const std::size_t n = model_.num_dofs();
const auto& elems = model_.elements();
std::vector<double> f(n, 0.0);
// Equivalent nodal load f_e = E*A*eps0 * b, b = (-c, -s, c, s). A member
// with eps0 > 0 (too long / heated) pushes its ends apart.
for (std::size_t e = 0; e < elems.size(); ++e) {
if (!active[e]) continue;
const Element& el = elems[e];
const double eps0 = el.prestrain + el.alpha * el.dT;
if (eps0 == 0.0) continue;
const Vec2 dir = model_.element_direction(e);
const double coef = el.section.E * el.section.A * eps0;
const auto g = element_dofs(el);
f[g[0]] += coef * (-dir.x);
f[g[1]] += coef * (-dir.y);
f[g[2]] += coef * (dir.x);
f[g[3]] += coef * (dir.y);
}
return f;
}
SolveResult StaticSolver::solve() const {
model_.validate();
const std::size_t n = model_.num_dofs();
const std::size_t ne = model_.num_elements();
const auto& elems = model_.elements();
const auto& springs = model_.springs();
// Precompute per-element geometry / installed strain.
std::vector<double> L0(ne), cc(ne), ss(ne), eps0(ne), EA(ne);
for (std::size_t e = 0; e < ne; ++e) {
const Element& el = elems[e];
L0[e] = model_.element_length(e);
const Vec2 dir = model_.element_direction(e);
cc[e] = dir.x;
ss[e] = dir.y;
eps0[e] = el.prestrain + el.alpha * el.dT;
EA[e] = el.section.E * el.section.A;
}
// Load-independent part of the right-hand side: applied loads plus the
// equivalent load from prescribed spring settlement (kappa * settlement *
// dir at the supported node).
std::vector<double> f_base(n, 0.0);
for (const Load& ld : model_.loads()) {
f_base[ld.dof] += ld.value;
}
for (const SpringSupport& sp : springs) {
if (sp.settlement != 0.0) {
f_base[2 * sp.node] += sp.kappa * sp.settlement * sp.dir.x;
f_base[2 * sp.node + 1] += sp.kappa * sp.settlement * sp.dir.y;
}
}
// Total axial force of element e implied by displacement field u.
auto axial_from_u = [&](std::size_t e, const std::vector<double>& u) {
const Element& el = elems[e];
const double d1x = u[2 * el.n1];
const double d1y = u[2 * el.n1 + 1];
const double d2x = u[2 * el.n2];
const double d2y = u[2 * el.n2 + 1];
const double delta = cc[e] * (d2x - d1x) + ss[e] * (d2y - d1y);
return EA[e] * (delta / L0[e] - eps0[e]);
};
// A unilateral member is judged by its total axial force: a cable
// participates only in tension, a strut only in compression, a bar always.
auto should_be_active = [&](std::size_t e, double N) {
switch (elems[e].kind) {
case MemberKind::Cable: return N > 0.0;
case MemberKind::Strut: return N < 0.0;
case MemberKind::Bar: return true;
}
return true;
};
std::vector<bool> active(ne, true);
std::vector<double> N(ne, 0.0);
std::vector<double> u(n, 0.0);
// Iterate to a state where the assumed member forces (driving the geometric
// stiffness), the displacements they produce, and the forces recovered from
// those displacements all coincide, with a consistent unilateral active set.
const int kMaxIter = 1000;
for (int iter = 0; iter < kMaxIter; ++iter) {
DenseMatrix K = assemble_system(active);
// Initial-stress (geometric) stiffness of active members using the
// current assumed force: N/L0 * b_t^T b_t, b_t = (s, -c, -s, c).
for (std::size_t e = 0; e < ne; ++e) {
if (!active[e] || N[e] == 0.0) continue;
const double factor = N[e] / L0[e];
const double bt[4] = {ss[e], -cc[e], -ss[e], cc[e]};
const auto g = element_dofs(elems[e]);
for (std::size_t i = 0; i < 4; ++i) {
for (std::size_t j = 0; j < 4; ++j) {
K(g[i], g[j]) += factor * bt[i] * bt[j];
}
}
}
std::vector<double> f = f_base;
const std::vector<double> fp = prestrain_force(active);
for (std::size_t i = 0; i < n; ++i) f[i] += fp[i];
const std::vector<double> u_new = solve_spd(K, f);
// Recover candidate forces for every member from the new displacements
// and decide the consistent active set.
std::vector<double> N_new(ne);
std::vector<bool> active_new(ne);
double scale = 0.0;
for (std::size_t e = 0; e < ne; ++e) {
N_new[e] = axial_from_u(e, u_new);
active_new[e] = should_be_active(e, N_new[e]);
scale = std::max(scale, std::fabs(N_new[e]));
}
const bool active_same = (active_new == active);
double max_dN = 0.0;
for (std::size_t e = 0; e < ne; ++e) {
if (active[e] && active_new[e]) {
max_dN = std::max(max_dN, std::fabs(N_new[e] - N[e]));
}
}
const bool force_converged = max_dN <= 1e-10 * scale + 1e-9;
u = u_new;
N = N_new;
active = active_new;
if (active_same && force_converged) break;
}
SolveResult result;
result.displacements = u;
result.reactions.assign(n, 0.0);
result.axial_forces.assign(ne, 0.0);
result.utilization.assign(ne, 0.0);
result.active = active;
for (std::size_t e = 0; e < ne; ++e) {
if (!active[e]) continue; // dropped member: exactly zero force / util.
result.axial_forces[e] = N[e];
result.utilization[e] = N[e] / (kAllowableStress * elems[e].section.A);
}
// Support reactions: the force the spring applies to the structure, which
// opposes the spring extension measured relative to any settlement offset.
for (const SpringSupport& sp : springs) {
const double along =
u[2 * sp.node] * sp.dir.x + u[2 * sp.node + 1] * sp.dir.y;
const double extension = along - sp.settlement;
const double force = -sp.kappa * extension;
result.reactions[2 * sp.node] += force * sp.dir.x;
result.reactions[2 * sp.node + 1] += force * sp.dir.y;
}
return result;
}
} // namespace truss2d
/opt/truss2d/src/solver.cpp
#include <array> #include <cmath>
#include <algorithm> #include <array> #include <cmath>
Configure, build, and run tests
cmake -S /opt/truss2d -B /opt/truss2d/build -DCMAKE_BUILD_TYPE=Release >/dev/null 2>&1 && cmake --build /opt/truss2d/build --parallel 2>&1 | tail -20 && cd /opt/truss2d/build && ctest --output-on-failure 2>&1 | tail -20
[ 11%] Building CXX object CMakeFiles/truss2d.dir/src/dense_matrix.cpp.o
[ 22%] Building CXX object CMakeFiles/truss2d.dir/src/model.cpp.o
[ 33%] Building CXX object CMakeFiles/truss2d.dir/src/io.cpp.o
[ 44%] Building CXX object CMakeFiles/truss2d.dir/src/solver.cpp.o
[ 55%] Linking CXX static library libtruss2d.a
[ 55%] Built target truss2d
[ 66%] Building CXX object CMakeFiles/truss2d_cli.dir/src/main.cpp.o
[ 77%] Building CXX object tests/CMakeFiles/test_smoke.dir/test_smoke.cpp.o
[ 88%] Linking CXX executable truss2d_cli
[ 88%] Built target truss2d_cli
[100%] Linking CXX executable test_smoke
[100%] Built target test_smoke
Test project /opt/truss2d/build
Start 1: smoke
1/1 Test #1: smoke ............................ Passed 0.00 sec
100% tests passed, 0 tests failed out of 1
Total Test time (real) = 0.00 sec
[stdout]
[ 11%] Building CXX object CMakeFiles/truss2d.dir/src/dense_matrix.cpp.o
[ 22%] Building CXX object CMakeFiles/truss2d.dir/src/model.cpp.o
[ 33%] Building CXX object CMakeFiles/truss2d.dir/src/io.cpp.o
[ 44%] Building CXX object CMakeFiles/truss2d.dir/src/solver.cpp.o
[ 55%] Linking CXX static library libtruss2d.a
[ 55%] Built target truss2d
[ 66%] Building CXX object CMakeFiles/truss2d_cli.dir/src/main.cpp.o
[ 77%] Building CXX object tests/CMakeFiles/test_smoke.dir/test_smoke.cpp.o
[ 88%] Linking CXX executable truss2d_cli
[ 88%] Built target truss2d_cli
[100%] Linking CXX executable test_smoke
[100%] Built target test_smoke
Test project /opt/truss2d/build
Start 1: smoke
1/1 Test #1: smoke ............................ Passed 0.00 sec
100% tests passed, 0 tests failed out of 1
Total Test time (real) = 0.00 secCompile and run analytic verification
cat > /tmp/verify.cpp <<'EOF'
#include <cmath>
#include <cstdio>
#include <vector>
#include "truss2d/model.hpp"
#include "truss2d/solver.hpp"
using namespace truss2d;
static int fails=0;
void chk(const char* name,double got,double want,double tol){
bool ok=std::fabs(got-want)<=tol+1e-9*std::fabs(want);
printf("[%s] %-30s got=%.6e want=%.6e %s\n",ok?"PASS":"FAIL",name,got,want,ok?"":" <<<");
if(!ok)++fails;
}
int main(){
// ---- Test 1: single horizontal bar, axial load, stiff transverse springs ----
// Node0 fully sprung (very stiff), node1 sprung in y. Pull node1 +x by P.
// bar EA/L stiffness; u1x = P/(EA/L) approx (springs huge so node0 ~fixed).
{
Model m; double E=2e11,A=1e-3,L=2.0; double EA=E*A;
m.add_node(0,0); m.add_node(L,0);
m.add_element(0,1,Section{E,A});
double K=1e15;
m.add_spring(0,1,0,K); m.add_spring(0,0,1,K); m.add_spring(1,0,1,K);
double P=1e4; m.add_load(2,P); // node1 x
auto r=StaticSolver(m).solve();
double u1x=r.displacements[2];
double expect_u=P/(EA/L); // node0 nearly fixed
chk("bar u1x",u1x,expect_u,1e-9);
chk("bar axial",r.axial_forces[0],P,1.0); // tension = P
chk("bar util",r.utilization[0],P/(2.5e8*A),1e-6);
// reaction at node0 x should be -P
chk("bar React0x",r.reactions[0],-P,1.0);
}
// ---- Test 2: prestrain compressive (too long) ----
// Single bar between two stiffly sprung nodes, prestrain eps>0 -> compression N=-EA*eps
{
Model m; double E=2e11,A=1e-3,L=2.0; double EA=E*A; double eps=1e-4;
m.add_node(0,0); m.add_node(L,0);
m.add_element(0,1,Section{E,A},MemberKind::Bar,eps);
double K=1e16;
m.add_spring(0,1,0,K);m.add_spring(0,0,1,K);m.add_spring(1,1,0,K);m.add_spring(1,0,1,K);
auto r=StaticSolver(m).solve();
chk("prestrain axial",r.axial_forces[0],-EA*eps,1.0); // compression
}
// ---- Test 3: thermal expansion restrained -> compression ----
{
Model m; double E=2e11,A=1e-3,L=2.0; double EA=E*A; double alpha=1.2e-5,dT=50;
m.add_node(0,0); m.add_node(L,0);
m.add_element(0,1,Section{E,A},MemberKind::Bar,0.0,alpha,dT);
double K=1e16;
m.add_spring(0,1,0,K);m.add_spring(0,0,1,K);m.add_spring(1,1,0,K);m.add_spring(1,0,1,K);
auto r=StaticSolver(m).solve();
chk("thermal axial",r.axial_forces[0],-EA*alpha*dT,1.0);
}
// ---- Test 4: cable goes slack under compression load ----
// horizontal cable; push node1 toward node0 (compressive) -> slack, N=0, inactive
{
Model m; double E=2e11,A=1e-3,L=2.0;
m.add_node(0,0); m.add_node(L,0);
m.add_element(0,1,Section{E,A},MemberKind::Cable);
double K=1e15;
m.add_spring(0,1,0,K);m.add_spring(0,0,1,K);m.add_spring(1,0,1,K);
m.add_load(2,-1e4); // push node1 in -x (toward node0) => compression
auto r=StaticSolver(m).solve();
chk("cable slack axial",r.axial_forces[0],0.0,1e-6);
printf(" cable active=%d (want 0)\n",(int)r.active[0]);
if(r.active[0])++fails;
}
// ---- Test 5: cable taut under tension ----
{
Model m; double E=2e11,A=1e-3,L=2.0; double EA=E*A;
m.add_node(0,0); m.add_node(L,0);
m.add_element(0,1,Section{E,A},MemberKind::Cable);
double K=1e15;
m.add_spring(0,1,0,K);m.add_spring(0,0,1,K);m.add_spring(1,0,1,K);
double P=1e4; m.add_load(2,P);
auto r=StaticSolver(m).solve();
chk("cable taut axial",r.axial_forces[0],P,1.0);
printf(" cable active=%d (want 1)\n",(int)r.active[0]);
if(!r.active[0])++fails;
}
// ---- Test 6: strut separates under tension load ----
{
Model m; double E=2e11,A=1e-3,L=2.0;
m.add_node(0,0); m.add_node(L,0);
m.add_element(0,1,Section{E,A},MemberKind::Strut);
double K=1e15;
m.add_spring(0,1,0,K);m.add_spring(0,0,1,K);m.add_spring(1,0,1,K);
m.add_load(2,1e4); // pull apart => tension => strut separates
auto r=StaticSolver(m).solve();
chk("strut sep axial",r.axial_forces[0],0.0,1e-6);
if(r.active[0])++fails;
}
// ---- Test 7: settlement loads structure ----
// single spring with settlement, no load. Node only sprung in x with settlement s.
// bar to a fixed node. spring at node1 x with settlement delta.
// Actually test reaction: node fixed by stiff springs except one spring with settlement.
{
Model m; double E=2e11,A=1e-3,L=2.0;
m.add_node(0,0); m.add_node(L,0);
m.add_element(0,1,Section{E,A});
double Kstiff=1e16, Ksoft=1e6, sett=0.01;
m.add_spring(0,1,0,Kstiff);m.add_spring(0,0,1,Kstiff);m.add_spring(1,0,1,Kstiff);
m.add_spring(1,1,0,Ksoft,sett); // soft spring x at node1 with settlement
auto r=StaticSolver(m).solve();
// node0 ~fixed. bar stiff EA/L=1e8. soft spring 1e6 pulls node1 toward sett.
// series: node1 displacement u: bar force = EA/L*u ; spring force= Ksoft*(u-sett) opposing
// equilibrium node1 x: -EA/L*u - Ksoft*(u-sett)=0 -> u*(EA/L+Ksoft)=Ksoft*sett
double kbar=E*A/L;
double u1=Ksoft*sett/(kbar+Ksoft);
chk("settle u1x",r.displacements[2],u1,1e-9);
// reaction soft spring = -Ksoft*(u1-sett)
double Rsoft=-Ksoft*(u1-sett);
// node1 x reaction only from soft spring (stiff is in y)
chk("settle React1x",r.reactions[2],Rsoft,1e-3);
}
// ---- Test 8: geometric stiffness of pretensioned cable (transverse) ----
// Horizontal cable node0-node1-? Use a taut string: two cables meeting at center node,
// pretensioned, transverse load -> transverse stiffness ~ 2*T/L * ... Compare to linear (no geom) which would be ~0 transverse from axial alone?
// Simpler: single horizontal bar with large prestrain tension, apply transverse load at node1,
// node1 free in y but held by soft y-spring; geometric adds N/L transverse stiffness.
{
Model m; double E=2e11,A=1e-3,L=2.0;
m.add_node(0,0); m.add_node(L,0);
// tension via prestrain negative (too short -> tension). eps0<0 => N=-EA*eps0>0
double eps=-1e-3; double EA=E*A; double T=-EA*eps; // tension
m.add_element(0,1,Section{E,A},MemberKind::Bar,eps);
double Kstiff=1e16;
m.add_spring(0,1,0,Kstiff);m.add_spring(0,0,1,Kstiff);m.add_spring(1,1,0,Kstiff);
double Ksofty=1.0; // tiny y spring at node1 so it's supported
m.add_spring(1,0,1,Ksofty);
double Pt=100.0; m.add_load(3,Pt); // transverse y at node1
auto r=StaticSolver(m).solve();
// node1 held in x (stiff), free-ish in y. Transverse stiffness = geom N/L + softy.
// equilibrium: (T/L + Ksofty)*uy = Pt approximately (geom from member tension T)
double uy=r.displacements[3];
double kgeo=T/L;
double expect_uy=Pt/(kgeo+Ksofty);
chk("geom transverse uy",uy,expect_uy,1e-6*std::fabs(expect_uy)+1e-9);
printf(" T=%.3e kgeo=%.3e uy=%.6e expect=%.6e\n",T,kgeo,uy,expect_uy);
}
printf("\n%s (%d failures)\n", fails==0?"ALL PASS":"SOME FAILED", fails);
return fails?1:0;
}
EOF
g++ -std=c++17 -I/opt/truss2d/include /tmp/verify.cpp /opt/truss2d/build/libtruss2d.a -o /tmp/verify 2>&1 | head -30 && /tmp/verifyExit code 134 terminate called after throwing an instance of 'std::runtime_error' what(): solve_spd: matrix is singular or not positive-definite /bin/bash: line 311: 1810 Aborted /tmp/verify [error] tool reported failure
#include "truss2d/solver.hpp"
#include <algorithm>
#include <array>
#include <cmath>
#include <cstddef>
#include <stdexcept>
#include <vector>
#include "truss2d/dense_matrix.hpp"
#include "truss2d/types.hpp"
namespace truss2d {
namespace {
// Per-element kinematic/material data precomputed once for solve().
//
// Local DOF order is {n1.x, n1.y, n2.x, n2.y}. With direction cosines (c, s)
// from n1 toward n2:
// b = [-c, -s, c, s] maps nodal displacements to axial extension: e = b . u
// g = [ s, -c, -s, c] maps them to the transverse relative displacement.
// b and g are orthogonal: b is the member axis, g the perpendicular, so the
// geometric (initial-stress) term built from g vanishes for axial motion and
// for rigid translation of both ends.
struct ElemData {
std::array<std::size_t, 4> dof;
double L0; // undeformed length
double EA; // axial rigidity
double eps0; // combined installed (mechanical + thermal) strain
std::array<double, 4> b;
std::array<double, 4> g;
};
// Allowable axial stress [Pa]; member capacity scales with area.
constexpr double kSigmaAllow = 2.5e8;
} // namespace
DenseMatrix StaticSolver::element_stiffness_global(std::size_t e) const {
// element_length / element_direction validate the index and zero length.
const double L0 = model_.element_length(e);
const Vec2 d = model_.element_direction(e);
const Element& el = model_.elements()[e];
const double k = el.section.E * el.section.A / L0;
const double c = d.x, s = d.y;
const double b[4] = {-c, -s, c, s};
DenseMatrix Ke(4, 4);
for (int i = 0; i < 4; ++i) {
for (int j = 0; j < 4; ++j) {
Ke(i, j) = k * b[i] * b[j];
}
}
return Ke;
}
DenseMatrix StaticSolver::assemble_system(const std::vector<bool>& active) const {
const std::size_t n = model_.num_dofs();
DenseMatrix K(n, n);
const auto& elems = model_.elements();
for (std::size_t e = 0; e < elems.size(); ++e) {
const bool act = (e >= active.size()) || active[e];
if (!act) continue;
const DenseMatrix Ke = element_stiffness_global(e);
const Element& el = elems[e];
const std::size_t dof[4] = {2 * el.n1, 2 * el.n1 + 1,
2 * el.n2, 2 * el.n2 + 1};
for (int i = 0; i < 4; ++i) {
for (int j = 0; j < 4; ++j) {
K(dof[i], dof[j]) += Ke(i, j);
}
}
}
// Grounded spring supports: each adds kappa * (dir outer dir) at its node's
// 2x2 block. A spring resists only the displacement component along `dir`.
for (const auto& sp : model_.springs()) {
const std::size_t dx = 2 * sp.node, dy = 2 * sp.node + 1;
const double kx = sp.kappa * sp.dir.x;
const double ky = sp.kappa * sp.dir.y;
K(dx, dx) += kx * sp.dir.x;
K(dx, dy) += kx * sp.dir.y;
K(dy, dx) += ky * sp.dir.x;
K(dy, dy) += ky * sp.dir.y;
}
return K;
}
std::vector<double>
StaticSolver::prestrain_force(const std::vector<bool>& active) const {
const std::size_t n = model_.num_dofs();
std::vector<double> F(n, 0.0);
const auto& elems = model_.elements();
for (std::size_t e = 0; e < elems.size(); ++e) {
const bool act = (e >= active.size()) || active[e];
if (!act) continue;
const double L0 = model_.element_length(e);
const Vec2 d = model_.element_direction(e);
const Element& el = elems[e];
const double EA = el.section.E * el.section.A;
const double eps0 = el.prestrain + el.alpha * el.dT;
const double c = d.x, s = d.y;
const double b[4] = {-c, -s, c, s};
const std::size_t dof[4] = {2 * el.n1, 2 * el.n1 + 1,
2 * el.n2, 2 * el.n2 + 1};
// Equivalent nodal force from installed strain: K u = F_ext + EA*eps0*b.
// eps0 > 0 (too long / heated) pushes the ends apart -> compression.
for (int i = 0; i < 4; ++i) {
F[dof[i]] += EA * eps0 * b[i];
}
}
return F;
}
SolveResult StaticSolver::solve() const {
model_.validate();
const std::size_t n = model_.num_dofs();
const auto& elems = model_.elements();
const std::size_t ne = elems.size();
// --- Precompute per-element data ---------------------------------------
std::vector<ElemData> ed(ne);
for (std::size_t e = 0; e < ne; ++e) {
const Element& el = elems[e];
const double L0 = model_.element_length(e);
const Vec2 d = model_.element_direction(e);
const double c = d.x, s = d.y;
ElemData x;
x.dof = {2 * el.n1, 2 * el.n1 + 1, 2 * el.n2, 2 * el.n2 + 1};
x.L0 = L0;
x.EA = el.section.E * el.section.A;
x.eps0 = el.prestrain + el.alpha * el.dT;
x.b = {-c, -s, c, s};
x.g = {s, -c, -s, c};
ed[e] = x;
}
// --- Constant right-hand-side contributions ----------------------------
std::vector<double> Fapplied(n, 0.0);
for (const auto& ld : model_.loads()) {
Fapplied[ld.dof] += ld.value;
}
// Support settlement: a prescribed grounded-end offset loads the structure
// even with no applied force. (K_spring) u = ... + kappa*settlement*dir.
std::vector<double> Fsettle(n, 0.0);
for (const auto& sp : model_.springs()) {
Fsettle[2 * sp.node] += sp.kappa * sp.settlement * sp.dir.x;
Fsettle[2 * sp.node + 1] += sp.kappa * sp.settlement * sp.dir.y;
}
// Total axial force of a member from a displacement field (tension +).
auto member_force = [&](std::size_t e, const std::vector<double>& u) {
const ElemData& x = ed[e];
double ext = 0.0;
for (int i = 0; i < 4; ++i) ext += x.b[i] * u[x.dof[i]];
return x.EA * (ext / x.L0 - x.eps0);
};
std::vector<bool> active(ne, true);
std::vector<double> u(n, 0.0);
std::vector<double> N(ne, 0.0);
std::vector<std::vector<bool>> seen; // active-set cycle guard
const int kMaxOuter = static_cast<int>(ne) + 64;
for (int outer = 0; outer < kMaxOuter; ++outer) {
// --- Inner loop: consistent equilibrium with initial-stress effect --
//
// The transverse (geometric) stiffness of a taut member depends on the
// force it carries, and that force depends on the displacement. We seek
// the state where the assumed forces, the displacements they produce,
// and the forces recovered from those displacements all coincide. The
// geometric term keeps the matrix symmetric, so a fixed-point (Picard)
// iteration paired with the supplied SPD solver is the natural scheme;
// it converges fast because the initial-stress term is a small
// correction to the elastic stiffness.
std::vector<double> Ncur(ne, 0.0);
for (int inner = 0; inner < 100; ++inner) {
DenseMatrix K = assemble_system(active);
for (std::size_t e = 0; e < ne; ++e) {
if (!active[e]) continue;
const double scale = Ncur[e] / ed[e].L0;
if (scale == 0.0) continue;
const ElemData& x = ed[e];
for (int i = 0; i < 4; ++i) {
for (int j = 0; j < 4; ++j) {
K(x.dof[i], x.dof[j]) += scale * x.g[i] * x.g[j];
}
}
}
std::vector<double> F = prestrain_force(active);
for (std::size_t i = 0; i < n; ++i) F[i] += Fapplied[i] + Fsettle[i];
u = solve_spd(K, F); // propagates on singular / indefinite systems
std::vector<double> Nnew(ne, 0.0);
double maxchange = 0.0, maxmag = 0.0;
for (std::size_t e = 0; e < ne; ++e) {
Nnew[e] = member_force(e, u);
if (active[e]) {
maxchange = std::max(maxchange, std::abs(Nnew[e] - Ncur[e]));
maxmag = std::max(maxmag, std::abs(Nnew[e]));
}
}
Ncur.swap(Nnew);
if (maxchange <= 1e-12 * maxmag + 1e-12) break;
}
N = Ncur;
// --- Active-set update for unilateral members ----------------------
// A cable carries tension only; a strut compression only. Each is
// judged by its total axial force (installed + elastic). An inactive
// member's force here is the value it would carry if reinstated.
bool changed = false;
std::vector<bool> newactive = active;
for (std::size_t e = 0; e < ne; ++e) {
const MemberKind kind = elems[e].kind;
if (kind == MemberKind::Bar) continue;
const bool want = (kind == MemberKind::Cable) ? (N[e] > 0.0)
: (N[e] < 0.0);
if (want != active[e]) {
changed = true;
newactive[e] = want;
}
}
if (!changed) break;
// Guard against a (theoretically excluded) active-set cycle.
if (std::find(seen.begin(), seen.end(), newactive) != seen.end()) break;
seen.push_back(active);
active = newactive;
}
// --- Assemble outputs --------------------------------------------------
SolveResult res;
res.displacements = u;
res.reactions.assign(n, 0.0);
for (const auto& sp : model_.springs()) {
const double along =
u[2 * sp.node] * sp.dir.x + u[2 * sp.node + 1] * sp.dir.y;
const double ext = along - sp.settlement;
const double fmag = sp.kappa * ext; // spring extension force
// Force the spring applies to the structure opposes its extension.
res.reactions[2 * sp.node] += -fmag * sp.dir.x;
res.reactions[2 * sp.node + 1] += -fmag * sp.dir.y;
}
res.axial_forces.assign(ne, 0.0);
res.utilization.assign(ne, 0.0);
res.active.assign(ne, true);
for (std::size_t e = 0; e < ne; ++e) {
const bool act =
(elems[e].kind == MemberKind::Bar) ? true : active[e];
res.active[e] = act;
if (act) {
res.axial_forces[e] = N[e];
res.utilization[e] = N[e] / (kSigmaAllow * elems[e].section.A);
}
// Dropped unilateral members report exactly zero force / utilization.
}
return res;
}
} // namespace truss2d
/opt/truss2d/src/solver.cpp
1 #include "truss2d/solver.hpp"
2
3 #include <algorithm>
4 #include <array>
5 #include <cmath>
6 #include <cstddef>
7 #include <stdexcept>
8 #include <vector>
9
10 #include "truss2d/dense_matrix.hpp"
11 #include "truss2d/types.hpp"
12
13 namespace truss2d {
14
15 // ============================================================================
16 // Static solver for a 2D network of two-force members (bars/cables/struts) on
17 // finite spring supports, with installed (mechanical + thermal) strain and the
18 // small-displacement initial-stress (geometric) stiffening of taut members.
19 //
20 // Kinematics (small displacement). For element e with direction cosines
21 // (c, s) from n1 -> n2 and undeformed length L0, the axial extension is the
22 // first-order change of end separation projected onto the member axis:
23 //
24 // delta = c*(u2x - u1x) + s*(u2y - u1y) = b . u_local,
25 // b = (-c, -s, c, s).
26 //
27 // Axial force (tension positive). The installed free strain is
28 // eps0 = prestrain + alpha*dT
29 // (a member that is too long / heated has eps0 > 0 and is compressive when
30 // restrained). The elastic strain is the actual extension minus the free
31 // strain, so
32 // N = E*A * (delta / L0 - eps0).
33 //
34 // Governing system. For a given active set,
35 // ( K_elastic(active) + K_spring + K_geom(N) ) u = f_load + f_settle
36 // + f_prestrain(active),
37 // where the geometric stiffness of an active member is N/L0 times the outer
38 // product of the transverse kinematic vector b_t = (s, -c, -s, c), so it acts
39 // only on relative transverse motion, vanishes for rigid translation and for
40 // motion along the axis, stiffens under tension and softens under compression.
41 // K_geom depends on N which depends on u, so solve() iterates the member forces
42 // to a self-consistent state. The matrix helpers below report the ordinary
43 // elastic behaviour only; the geometric term lives entirely inside solve().
44 // ============================================================================
45
46 namespace {
47
48 constexpr double kAllowableStress = 2.5e8; // [Pa]
49
50 // Global DOF indices for an element's local DOFs {n1.x, n1.y, n2.x, n2.y}.
51 std::array<std::size_t, 4> element_dofs(const Element& el) {
52 return {2 * el.n1, 2 * el.n1 + 1, 2 * el.n2, 2 * el.n2 + 1};
53 }
54
55 } // namespace
56
57 DenseMatrix StaticSolver::element_stiffness_global(std::size_t e) const {
58 const Element& el = model_.elements()[e];
59 const double L = model_.element_length(e);
60 const Vec2 dir = model_.element_direction(e);
61 const double c = dir.x;
62 const double s = dir.y;
63 const double k = el.section.E * el.section.A / L;
64
65 // K_e = (E A / L) * b^T b, with b = (-c, -s, c, s).
66 const double b[4] = {-c, -s, c, s};
67 DenseMatrix Ke(4, 4);
68 for (std::size_t i = 0; i < 4; ++i) {
69 for (std::size_t j = 0; j < 4; ++j) {
70 Ke(i, j) = k * b[i] * b[j];
71 }
72 }
73 return Ke;
74 }
75
76 DenseMatrix StaticSolver::assemble_system(const std::vector<bool>& active) const {
77 const std::size_t n = model_.num_dofs();
78 const auto& elems = model_.elements();
79 DenseMatrix K(n, n);
80
81 // Elastic stiffness of every participating member.
82 for (std::size_t e = 0; e < elems.size(); ++e) {
83 if (!active[e]) continue;
84 const DenseMatrix Ke = element_stiffness_global(e);
85 const auto g = element_dofs(elems[e]);
86 for (std::size_t i = 0; i < 4; ++i) {
87 for (std::size_t j = 0; j < 4; ++j) {
88 K(g[i], g[j]) += Ke(i, j);
89 }
90 }
91 }
92
93 // Spring supports: kappa * dir dir^T in the supported node's 2x2 block.
94 for (const SpringSupport& sp : model_.springs()) {
95 const std::size_t gx = 2 * sp.node;
96 const std::size_t gy = 2 * sp.node + 1;
97 const double kx = sp.kappa * sp.dir.x;
98 const double ky = sp.kappa * sp.dir.y;
99 K(gx, gx) += kx * sp.dir.x;
100 K(gx, gy) += kx * sp.dir.y;
101 K(gy, gx) += ky * sp.dir.x;
102 K(gy, gy) += ky * sp.dir.y;
103 }
104
105 return K;
106 }
107
108 std::vector<double>
109 StaticSolver::prestrain_force(const std::vector<bool>& active) const {
110 const std::size_t n = model_.num_dofs();
111 const auto& elems = model_.elements();
112 std::vector<double> f(n, 0.0);
113
114 // Equivalent nodal load f_e = E*A*eps0 * b, b = (-c, -s, c, s). A member
115 // with eps0 > 0 (too long / heated) pushes its ends apart.
116 for (std::size_t e = 0; e < elems.size(); ++e) {
117 if (!active[e]) continue;
118 const Element& el = elems[e];
119 const double eps0 = el.prestrain + el.alpha * el.dT;
120 if (eps0 == 0.0) continue;
121 const Vec2 dir = model_.element_direction(e);
122 const double coef = el.section.E * el.section.A * eps0;
123 const auto g = element_dofs(el);
124 f[g[0]] += coef * (-dir.x);
125 f[g[1]] += coef * (-dir.y);
126 f[g[2]] += coef * (dir.x);
127 f[g[3]] += coef * (dir.y);
128 }
129 return f;
130 }
131
132 SolveResult StaticSolver::solve() const {
133 model_.validate();
134
135 const std::size_t n = model_.num_dofs();
136 const std::size_t ne = model_.num_elements();
137 const auto& elems = model_.elements();
138 const auto& springs = model_.springs();
139
140 // Precompute per-element geometry / installed strain.
141 std::vector<double> L0(ne), cc(ne), ss(ne), eps0(ne), EA(ne);
142 for (std::size_t e = 0; e < ne; ++e) {
143 const Element& el = elems[e];
144 L0[e] = model_.element_length(e);
145 const Vec2 dir = model_.element_direction(e);
146 cc[e] = dir.x;
147 ss[e] = dir.y;
148 eps0[e] = el.prestrain + el.alpha * el.dT;
149 EA[e] = el.section.E * el.section.A;
150 }
151
152 // Load-independent part of the right-hand side: applied loads plus the
153 // equivalent load from prescribed spring settlement (kappa * settlement *
154 // dir at the supported node).
155 std::vector<double> f_base(n, 0.0);
156 for (const Load& ld : model_.loads()) {
157 f_base[ld.dof] += ld.value;
158 }
159 for (const SpringSupport& sp : springs) {
160 if (sp.settlement != 0.0) {
161 f_base[2 * sp.node] += sp.kappa * sp.settlement * sp.dir.x;
162 f_base[2 * sp.node + 1] += sp.kappa * sp.settlement * sp.dir.y;
163 }
164 }
165
166 // Total axial force of element e implied by displacement field u.
167 auto axial_from_u = [&](std::size_t e, const std::vector<double>& u) {
168 const Element& el = elems[e];
169 const double d1x = u[2 * el.n1];
170 const double d1y = u[2 * el.n1 + 1];
171 const double d2x = u[2 * el.n2];
172 const double d2y = u[2 * el.n2 + 1];
173 const double delta = cc[e] * (d2x - d1x) + ss[e] * (d2y - d1y);
174 return EA[e] * (delta / L0[e] - eps0[e]);
175 };
176
177 // A unilateral member is judged by its total axial force: a cable
178 // participates only in tension, a strut only in compression, a bar always.
179 auto should_be_active = [&](std::size_t e, double N) {
180 switch (elems[e].kind) {
181 case MemberKind::Cable: return N > 0.0;
182 case MemberKind::Strut: return N < 0.0;
183 case MemberKind::Bar: return true;
184 }
185 return true;
186 };
187
188 std::vector<bool> active(ne, true);
189 std::vector<double> N(ne, 0.0);
190 std::vector<double> u(n, 0.0);
191
192 // Iterate to a state where the assumed member forces (driving the geometric
193 // stiffness), the displacements they produce, and the forces recovered from
194 // those displacements all coincide, with a consistent unilateral active set.
195 const int kMaxIter = 1000;
196 for (int iter = 0; iter < kMaxIter; ++iter) {
197 DenseMatrix K = assemble_system(active);
198
199 // Initial-stress (geometric) stiffness of active members using the
200 // current assumed force: N/L0 * b_t^T b_t, b_t = (s, -c, -s, c).
201 for (std::size_t e = 0; e < ne; ++e) {
202 if (!active[e] || N[e] == 0.0) continue;
203 const double factor = N[e] / L0[e];
204 const double bt[4] = {ss[e], -cc[e], -ss[e], cc[e]};
205 const auto g = element_dofs(elems[e]);
206 for (std::size_t i = 0; i < 4; ++i) {
207 for (std::size_t j = 0; j < 4; ++j) {
208 K(g[i], g[j]) += factor * bt[i] * bt[j];
209 }
210 }
211 }
212
213 std::vector<double> f = f_base;
214 const std::vector<double> fp = prestrain_force(active);
215 for (std::size_t i = 0; i < n; ++i) f[i] += fp[i];
216
217 const std::vector<double> u_new = solve_spd(K, f);
218
219 // Recover candidate forces for every member from the new displacements
220 // and decide the consistent active set.
221 std::vector<double> N_new(ne);
222 std::vector<bool> active_new(ne);
223 double scale = 0.0;
224 for (std::size_t e = 0; e < ne; ++e) {
225 N_new[e] = axial_from_u(e, u_new);
226 active_new[e] = should_be_active(e, N_new[e]);
227 scale = std::max(scale, std::fabs(N_new[e]));
228 }
229
230 const bool active_same = (active_new == active);
231 double max_dN = 0.0;
232 for (std::size_t e = 0; e < ne; ++e) {
233 if (active[e] && active_new[e]) {
234 max_dN = std::max(max_dN, std::fabs(N_new[e] - N[e]));
235 }
236 }
237 const bool force_converged = max_dN <= 1e-10 * scale + 1e-9;
238
239 u = u_new;
240 N = N_new;
241 active = active_new;
242
243 if (active_same && force_converged) break;
244 }
245
246 SolveResult result;
247 result.displacements = u;
248 result.reactions.assign(n, 0.0);
249 result.axial_forces.assign(ne, 0.0);
250 result.utilization.assign(ne, 0.0);
251 result.active = active;
252
253 for (std::size_t e = 0; e < ne; ++e) {
254 if (!active[e]) continue; // dropped member: exactly zero force / util.
255 result.axial_forces[e] = N[e];
256 result.utilization[e] = N[e] / (kAllowableStress * elems[e].section.A);
257 }
258
259 // Support reactions: the force the spring applies to the structure, which
260 // opposes the spring extension measured relative to any settlement offset.
261 for (const SpringSupport& sp : springs) {
262 const double along =
263 u[2 * sp.node] * sp.dir.x + u[2 * sp.node + 1] * sp.dir.y;
264 const double extension = along - sp.settlement;
265 const double force = -sp.kappa * extension;
266 result.reactions[2 * sp.node] += force * sp.dir.x;
267 result.reactions[2 * sp.node + 1] += force * sp.dir.y;
268 }
269
270 return result;
271 }
272
273 } // namespace truss2d
274
Patch test supports
cd /tmp
# Add soft backup x-spring at node1 in tests 4,5,6 so structure stays well-posed when member drops
sed -i 's| m.add_load(2,-1e4); // push node1 in -x (toward node0) => compression| m.add_spring(1,1,0,1.0); // soft backup\n m.add_load(2,-1e4);|' verify.cpp
sed -i 's| double P=1e4; m.add_load(2,P);\n auto r=StaticSolver(m).solve();\n chk("cable taut|XX|' verify.cpp
echo donedone Shell cwd was reset to /opt/truss2d [stdout] done [stderr] Shell cwd was reset to /opt/truss2d
cmake -S . -B build -DCMAKE_BUILD_TYPE=Release >/tmp/cmake.log 2>&1 && cmake --build build --parallel >/tmp/build.log 2>&1 && echo "BUILD OK" && cd build && ctest --output-on-failure 2>&1 | tail -20
BUILD OK
Test project /opt/truss2d/build
Start 1: smoke
1/1 Test #1: smoke ............................ Passed 0.00 sec
100% tests passed, 0 tests failed out of 1
Total Test time (real) = 0.00 sec
[stdout]
BUILD OK
Test project /opt/truss2d/build
Start 1: smoke
1/1 Test #1: smoke ............................ Passed 0.00 sec
100% tests passed, 0 tests failed out of 1
Total Test time (real) = 0.00 secCompile and run full verification
cat > /tmp/verify.cpp <<'EOF'
#include <cmath>
#include <cstdio>
#include <vector>
#include "truss2d/model.hpp"
#include "truss2d/solver.hpp"
using namespace truss2d;
static int fails=0;
void chk(const char* name,double got,double want,double tol){
bool ok=std::fabs(got-want)<=tol+1e-7*std::fabs(want);
printf("[%s] %-22s got=%.6e want=%.6e\n",ok?"PASS":"FAIL",name,got,want);
if(!ok)++fails;
}
int main(){
try{ // T1 axial bar
Model m; double E=2e11,A=1e-3,L=2.0,EA=E*A;
m.add_node(0,0);m.add_node(L,0);m.add_element(0,1,Section{E,A});
double K=1e15; m.add_spring(0,1,0,K);m.add_spring(0,0,1,K);m.add_spring(1,0,1,K);
double P=1e4;m.add_load(2,P);
auto r=StaticSolver(m).solve();
chk("bar u1x",r.displacements[2],P/(EA/L),1e-9);
chk("bar axial",r.axial_forces[0],P,1.0);
chk("bar util",r.utilization[0],P/(2.5e8*A),1e-6);
chk("bar R0x",r.reactions[0],-P,1.0);
}catch(std::exception&e){printf("T1 throw %s\n",e.what());++fails;}
try{ // T2 prestrain compressive (too long, eps>0)
Model m;double E=2e11,A=1e-3,EA=E*A,eps=1e-4;
m.add_node(0,0);m.add_node(2,0);m.add_element(0,1,Section{E,A},MemberKind::Bar,eps);
double K=1e16;m.add_spring(0,1,0,K);m.add_spring(0,0,1,K);m.add_spring(1,1,0,K);m.add_spring(1,0,1,K);
auto r=StaticSolver(m).solve();
chk("prestrain axial",r.axial_forces[0],-EA*eps,1.0);
}catch(std::exception&e){printf("T2 throw %s\n",e.what());++fails;}
try{ // T3 thermal expansion restrained -> compression
Model m;double E=2e11,A=1e-3,EA=E*A,al=1.2e-5,dT=50;
m.add_node(0,0);m.add_node(2,0);m.add_element(0,1,Section{E,A},MemberKind::Bar,0,al,dT);
double K=1e16;m.add_spring(0,1,0,K);m.add_spring(0,0,1,K);m.add_spring(1,1,0,K);m.add_spring(1,0,1,K);
auto r=StaticSolver(m).solve();
chk("thermal axial",r.axial_forces[0],-EA*al*dT,1.0);
}catch(std::exception&e){printf("T3 throw %s\n",e.what());++fails;}
try{ // T4 cable slack under compression
Model m;double E=2e11,A=1e-3;
m.add_node(0,0);m.add_node(2,0);m.add_element(0,1,Section{E,A},MemberKind::Cable);
double K=1e15;m.add_spring(0,1,0,K);m.add_spring(0,0,1,K);m.add_spring(1,0,1,K);
m.add_spring(1,1,0,1.0); // soft backup so well-posed when slack
m.add_load(2,-1e4);
auto r=StaticSolver(m).solve();
chk("cable slack N",r.axial_forces[0],0.0,1e-6);
if(r.active[0]){printf("FAIL cable should be inactive\n");++fails;}
}catch(std::exception&e){printf("T4 throw %s\n",e.what());++fails;}
try{ // T5 cable taut tension
Model m;double E=2e11,A=1e-3;
m.add_node(0,0);m.add_node(2,0);m.add_element(0,1,Section{E,A},MemberKind::Cable);
double K=1e15;m.add_spring(0,1,0,K);m.add_spring(0,0,1,K);m.add_spring(1,0,1,K);
m.add_spring(1,1,0,1.0);
m.add_load(2,1e4);
auto r=StaticSolver(m).solve();
chk("cable taut N",r.axial_forces[0],1e4,1.0);
if(!r.active[0]){printf("FAIL cable should be active\n");++fails;}
}catch(std::exception&e){printf("T5 throw %s\n",e.what());++fails;}
try{ // T6 strut separates under tension
Model m;double E=2e11,A=1e-3;
m.add_node(0,0);m.add_node(2,0);m.add_element(0,1,Section{E,A},MemberKind::Strut);
double K=1e15;m.add_spring(0,1,0,K);m.add_spring(0,0,1,K);m.add_spring(1,0,1,K);
m.add_spring(1,1,0,1.0);
m.add_load(2,1e4);
auto r=StaticSolver(m).solve();
chk("strut sep N",r.axial_forces[0],0.0,1e-6);
if(r.active[0]){printf("FAIL strut should be inactive\n");++fails;}
}catch(std::exception&e){printf("T6 throw %s\n",e.what());++fails;}
try{ // T7 settlement
Model m;double E=2e11,A=1e-3,L=2.0;
m.add_node(0,0);m.add_node(L,0);m.add_element(0,1,Section{E,A});
double Ks=1e16,Kso=1e6,sett=0.01;
m.add_spring(0,1,0,Ks);m.add_spring(0,0,1,Ks);m.add_spring(1,0,1,Ks);
m.add_spring(1,1,0,Kso,sett);
auto r=StaticSolver(m).solve();
double kbar=E*A/L,u1=Kso*sett/(kbar+Kso);
chk("settle u1x",r.displacements[2],u1,1e-9);
chk("settle R1x",r.reactions[2],-Kso*(u1-sett),1e-2);
}catch(std::exception&e){printf("T7 throw %s\n",e.what());++fails;}
try{ // T8 geometric transverse stiffening of pretensioned member
Model m;double E=2e11,A=1e-3,L=2.0,EA=E*A,eps=-1e-3,T=-EA*eps;
m.add_node(0,0);m.add_node(L,0);m.add_element(0,1,Section{E,A},MemberKind::Bar,eps);
double Ks=1e16;m.add_spring(0,1,0,Ks);m.add_spring(0,0,1,Ks);m.add_spring(1,1,0,Ks);
double Ksy=1.0;m.add_spring(1,0,1,Ksy);
double Pt=100.0;m.add_load(3,Pt);
auto r=StaticSolver(m).solve();
double uy=r.displacements[3],kgeo=T/L,expect=Pt/(kgeo+Ksy);
chk("geom transv uy",uy,expect,1e-4*std::fabs(expect));
printf(" T=%.3e kgeo=%.3e\n",T,kgeo);
}catch(std::exception&e){printf("T8 throw %s\n",e.what());++fails;}
try{ // T9 skew roller: 45-degree spring, check no perpendicular reaction
Model m;double E=2e11,A=1e-3;
m.add_node(0,0);m.add_node(3,4); // length5 diagonal
m.add_element(0,1,Section{E,A});
double K=1e16;
m.add_spring(0,1,0,K);m.add_spring(0,0,1,K); // node0 fixed
m.add_spring(1,0.6,0.8,1e7); // skew roller along member dir at node1
m.add_spring(1,-0.8,0.6,1e16); // strong perpendicular hold to be well posed
m.add_load(2,5e3);m.add_load(3,0);
auto r=StaticSolver(m).solve();
// just check it solves & reaction along skew opposes; sanity that reaction vector
// at node1 is finite
printf(" T9 R1=(%.3e,%.3e) u1=(%.3e,%.3e)\n",r.reactions[2],r.reactions[3],r.displacements[2],r.displacements[3]);
}catch(std::exception&e){printf("T9 throw %s\n",e.what());++fails;}
printf("\n%s (%d failures)\n",fails==0?"ALL PASS":"FAILED",fails);
return fails?1:0;
}
EOF
g++ -std=c++17 -I/opt/truss2d/include /tmp/verify.cpp /opt/truss2d/build/libtruss2d.a -o /tmp/verify 2>&1|head && /tmp/verify[PASS] bar u1x got=1.000000e-04 want=1.000000e-04 [PASS] bar axial got=1.000000e+04 want=1.000000e+04 [PASS] bar util got=4.000000e-02 want=4.000000e-02 [PASS] bar R0x got=-1.000000e+04 want=-1.000000e+04 [PASS] prestrain axial got=-2.000000e+04 want=-2.000000e+04 [PASS] thermal axial got=-1.200000e+05 want=-1.200000e+05 [PASS] cable slack N got=0.000000e+00 want=0.000000e+00 [PASS] cable taut N got=1.000000e+04 want=1.000000e+04 [PASS] strut sep N got=0.000000e+00 want=0.000000e+00 [PASS] settle u1x got=9.900990e-05 want=9.900990e-05 [PASS] settle R1x got=9.900990e+03 want=9.900990e+03 [PASS] geom transv uy got=9.999900e-04 want=9.999900e-04 T=2.000e+05 kgeo=1.000e+05 T9 R1=(-3.560e+03,1.920e+03) u1=(3.600e-05,4.800e-05) ALL PASS (0 failures) [stdout] [PASS] bar u1x got=1.000000e-04 want=1.000000e-04 [PASS] bar axial got=1.000000e+04 want=1.000000e+04 [PASS] bar util got=4.000000e-02 want=4.000000e-02 [PASS] bar R0x got=-1.000000e+04 want=-1.000000e+04 [PASS] prestrain axial got=-2.000000e+04 want=-2.000000e+04 [PASS] thermal axial got=-1.200000e+05 want=-1.200000e+05 [PASS] cable slack N got=0.000000e+00 want=0.000000e+00 [PASS] cable taut N got=1.000000e+04 want=1.000000e+04 [PASS] strut sep N got=0.000000e+00 want=0.000000e+00 [PASS] settle u1x got=9.900990e-05 want=9.900990e-05 [PASS] settle R1x got=9.900990e+03 want=9.900990e+03 [PASS] geom transv uy got=9.999900e-04 want=9.999900e-04 T=2.000e+05 kgeo=1.000e+05 T9 R1=(-3.560e+03,1.920e+03) u1=(3.600e-05,4.800e-05) ALL PASS (0 failures)
Full nodal equilibrium residual check
cat > /tmp/equil.cpp <<'EOF'
#include <cmath>
#include <cstdio>
#include <vector>
#include "truss2d/model.hpp"
#include "truss2d/solver.hpp"
using namespace truss2d;
// Reconstruct full nodal equilibrium: for each element, internal force on nodes
// = axial part + geometric transverse part, computed independently from u & N.
int main(){
Model m; double E=2e11;
// a small skewed network: 3 nodes, 2 cables + 1 bar, prestrain, a load, springs
m.add_node(0,0); m.add_node(4,0); m.add_node(2,3);
m.add_element(0,2,Section{E,1e-3},MemberKind::Bar,-5e-4); // pretensioned
m.add_element(1,2,Section{E,1.2e-3},MemberKind::Cable,-2e-4);
m.add_element(0,1,Section{E,8e-4},MemberKind::Bar,1e-4,1.2e-5,30);
m.add_spring(0,1,0,1e16); m.add_spring(0,0,1,1e16);
m.add_spring(1,0,1,1e16); m.add_spring(1,1,0,5e7,0.002); // skew/settle
m.add_load(4,-3e4); m.add_load(5,-2e4);
auto r=StaticSolver(m).solve();
std::size_t n=m.num_dofs();
std::vector<double> resid(n,0.0);
// applied loads
for(auto&ld:m.loads()) resid[ld.dof]+=ld.value;
// reactions (already force applied by springs to structure)
for(std::size_t i=0;i<n;++i) resid[i]+=r.reactions[i];
// member internal forces (force member exerts ON nodes) = -(axial b) - geom
for(std::size_t e=0;e<m.num_elements();++e){
if(!r.active[e])continue;
auto el=m.elements()[e];
double L=m.element_length(e); auto d=m.element_direction(e);
double c=d.x,s=d.y,N=r.axial_forces[e];
std::size_t g[4]={2*el.n1,2*el.n1+1,2*el.n2,2*el.n2+1};
double u[4]={r.displacements[g[0]],r.displacements[g[1]],r.displacements[g[2]],r.displacements[g[3]]};
// axial internal force on nodes: member in tension N pulls n1 toward n2 (+axis on n1)
double ax[4]={ c*N, s*N, -c*N, -s*N };
// geometric transverse: force = (N/L) bt (bt.u), bt=(s,-c,-s,c); this resists, exerted on nodes = -(N/L)bt(bt.u)
double bt[4]={s,-c,-s,c}; double w=0; for(int i=0;i<4;++i)w+=bt[i]*u[i];
for(int i=0;i<4;++i){ double fon = ax[i] - (N/L)*bt[i]*w; resid[g[i]]+=fon; }
}
double mx=0; for(double v:resid) mx=std::max(mx,std::fabs(v));
printf("max nodal equilibrium residual = %.4e\n",mx);
for(std::size_t e=0;e<m.num_elements();++e)
printf("elem %zu active=%d N=%.5e util=%.4e\n",e,(int)r.active[e],r.axial_forces[e],r.utilization[e]);
printf("%s\n", mx<1e-3?"EQUILIBRIUM OK":"EQUILIBRIUM FAIL");
return mx<1e-3?0:1;
}
EOF
g++ -std=c++17 -I/opt/truss2d/include /tmp/equil.cpp /opt/truss2d/build/libtruss2d.a -o /tmp/equil 2>&1|head && /tmp/equilmax nodal equilibrium residual = 3.8264e-11 elem 0 active=1 N=-3.90412e+04 util=-1.5616e-01 elem 1 active=1 N=1.50099e+04 util=5.0033e-02 elem 2 active=1 N=-1.50441e+02 util=-7.5221e-04 EQUILIBRIUM OK [stdout] max nodal equilibrium residual = 3.8264e-11 elem 0 active=1 N=-3.90412e+04 util=-1.5616e-01 elem 1 active=1 N=1.50099e+04 util=5.0033e-02 elem 2 active=1 N=-1.50441e+02 util=-7.5221e-04 EQUILIBRIUM OK
#include <cmath>
#include <cstdio>
#include <string>
#include "truss2d/model.hpp"
#include "truss2d/solver.hpp"
using namespace truss2d;
static int failures = 0;
static void chk(const std::string& name, double got, double want, double tol) {
double err = std::fabs(got - want);
double rel = err / (std::fabs(want) + 1e-30);
bool ok = err <= tol || rel <= 1e-6;
if (!ok) ++failures;
std::printf("%-40s got=% .8e want=% .8e %s\n", name.c_str(), got, want,
ok ? "OK" : " <<< FAIL");
}
int main() {
const double EA = 2.0e8; // E*A with E=2e11, A=1e-3
const double E = 2.0e11, A = 1.0e-3;
// ---- Case 1: axial bar with a finite x-spring at node0 ----------------
{
double L = 2.0, kappa0 = 1.0e6, P = 1.0e3;
Model m;
m.add_node(0, 0);
m.add_node(L, 0);
m.add_element(0, 1, Section{E, A}, MemberKind::Bar);
m.add_spring(0, 1, 0, kappa0);
m.add_spring(0, 0, 1, 1e14);
m.add_spring(1, 0, 1, 1e14);
m.add_load(2, P); // node1 x
auto r = StaticSolver(m).solve();
chk("axial: N", r.axial_forces[0], P, 1e-3);
chk("axial: u0x", r.displacements[0], P / kappa0, 1e-12);
chk("axial: u1x", r.displacements[2], P / kappa0 + P * L / EA, 1e-12);
chk("axial: reaction0x", r.reactions[0], -P, 1e-3);
chk("axial: util", r.utilization[0], P / (2.5e8 * A), 1e-9);
}
// ---- Case 2: prestrain (too-long member => compression) ---------------
{
double L = 2.0, kappa = 1.0e8, eps0 = 1.0e-4;
double a = EA * eps0 / (2 * EA / L + kappa);
double Nexp = EA * (2 * a / L - eps0);
Model m;
m.add_node(0, 0);
m.add_node(L, 0);
m.add_element(0, 1, Section{E, A}, MemberKind::Bar, eps0);
m.add_spring(0, 1, 0, kappa);
m.add_spring(1, 1, 0, kappa);
m.add_spring(0, 0, 1, 1e14);
m.add_spring(1, 0, 1, 1e14);
auto r = StaticSolver(m).solve();
chk("prestrain: N (compression<0)", r.axial_forces[0], Nexp, 1e-4);
chk("prestrain: u1x (+x)", r.displacements[2], a, 1e-12);
chk("prestrain: u0x (-x)", r.displacements[0], -a, 1e-12);
}
// ---- Case 2b: thermal equivalent (alpha*dT == prestrain) --------------
{
double L = 2.0, kappa = 1.0e8, alpha = 1.0e-5, dT = 10.0; // eps0=1e-4
double eps0 = alpha * dT;
double a = EA * eps0 / (2 * EA / L + kappa);
double Nexp = EA * (2 * a / L - eps0);
Model m;
m.add_node(0, 0);
m.add_node(L, 0);
m.add_element(0, 1, Section{E, A}, MemberKind::Bar, 0.0, alpha, dT);
m.add_spring(0, 1, 0, kappa);
m.add_spring(1, 1, 0, kappa);
m.add_spring(0, 0, 1, 1e14);
m.add_spring(1, 0, 1, 1e14);
auto r = StaticSolver(m).solve();
chk("thermal: N", r.axial_forces[0], Nexp, 1e-4);
}
// ---- Case 3: settlement loads the structure ---------------------------
{
// node1 x held by spring with prescribed settlement s0; node0 fixed.
double L = 2.0, kappa = 1.0e8, s0 = 1.0e-3;
Model m;
m.add_node(0, 0);
m.add_node(L, 0);
m.add_element(0, 1, Section{E, A}, MemberKind::Bar);
m.add_spring(0, 1, 0, 1e14);
m.add_spring(0, 0, 1, 1e14);
m.add_spring(1, 0, 1, 1e14);
m.add_spring(1, 1, 0, kappa, s0); // settling spring along x at node1
auto r = StaticSolver(m).solve();
// x-system: node0 fixed; node1: (EA/L + kappa) u1 = kappa*s0
double u1 = kappa * s0 / (EA / L + kappa);
double Nexp = EA * (u1 / L); // tension if pulled +x
chk("settle: u1x", r.displacements[2], u1, 1e-12);
chk("settle: N", r.axial_forces[0], Nexp, 1e-3);
// reaction at settling spring = -kappa*(u1 - s0)
chk("settle: reaction1x", r.reactions[2], -kappa * (u1 - s0), 1e-3);
}
// ---- Case 4: cable goes slack under compression -----------------------
{
// Two cables to a loaded node; push node so one cable would compress.
// Simple: single cable that would be compressed by a load -> drops.
double L = 2.0;
Model m;
m.add_node(0, 0);
m.add_node(L, 0);
m.add_element(0, 1, Section{E, A}, MemberKind::Cable);
m.add_spring(0, 1, 0, 1e14);
m.add_spring(0, 0, 1, 1e14);
m.add_spring(1, 1, 0, 1e6); // node1 held in x by a spring
m.add_spring(1, 0, 1, 1e14);
m.add_load(2, -1.0e3); // push node1 toward node0 -> cable compresses
auto r = StaticSolver(m).solve();
chk("cable slack: active", r.active[0] ? 1.0 : 0.0, 0.0, 0);
chk("cable slack: force=0", r.axial_forces[0], 0.0, 0);
chk("cable slack: util=0", r.utilization[0], 0.0, 0);
// With the cable gone, node1 x is carried only by its spring 1e6:
chk("cable slack: u1x", r.displacements[2], -1.0e3 / 1.0e6, 1e-12);
}
// ---- Case 4b: cable stays taut under tension --------------------------
{
double L = 2.0, P = 1.0e3;
Model m;
m.add_node(0, 0);
m.add_node(L, 0);
m.add_element(0, 1, Section{E, A}, MemberKind::Cable);
m.add_spring(0, 1, 0, 1e14);
m.add_spring(0, 0, 1, 1e14);
m.add_spring(1, 0, 1, 1e14);
m.add_load(2, P); // pull node1 +x -> tension
auto r = StaticSolver(m).solve();
chk("cable taut: active", r.active[0] ? 1.0 : 0.0, 1.0, 0);
chk("cable taut: N", r.axial_forces[0], P, 1e-3);
}
// ---- Case 5: geometric (initial-stress) transverse stiffening ---------
{
// node0 fixed; member pretensioned by prestrain eps0<0; node1 x heavily
// anchored so N ~= -EA*eps0; transverse stiffness at node1.y is
// (soft y-spring) + N/L.
double L = 2.0, eps0 = -1.0e-3, ks = 1.0e4, P = 1.0e3;
Model m;
m.add_node(0, 0);
m.add_node(L, 0);
m.add_element(0, 1, Section{E, A}, MemberKind::Bar, eps0);
m.add_spring(0, 1, 0, 1e14);
m.add_spring(0, 0, 1, 1e14);
m.add_spring(1, 1, 0, 1e14); // anchor node1 x -> sets tension
m.add_spring(1, 0, 1, ks); // soft transverse spring at node1 y
m.add_load(3, P); // node1 y load
auto r = StaticSolver(m).solve();
double N = r.axial_forces[0];
double u1y_exp = P / (ks + N / L); // geometric stiffening included
chk("geom: N (tension)", N, -EA * eps0, 5.0);
chk("geom: u1y w/ N/L stiffening", r.displacements[3], u1y_exp, 1e-9);
// Sanity: without stiffening it would be P/ks = 0.1, clearly different.
}
std::printf("\n%s (%d failures)\n", failures ? "SOME TESTS FAILED" : "ALL PASS",
failures);
return failures ? 1 : 0;
}
cd /opt/truss2d && g++ -std=c++17 -Iinclude verify.cpp src/model.cpp src/dense_matrix.cpp src/solver.cpp -o /tmp/verify 2>&1 | head -20 && /tmp/verify
axial: N got= 1.00000000e+03 want= 1.00000000e+03 OK axial: u0x got= 1.00000000e-03 want= 1.00000000e-03 OK axial: u1x got= 1.01000000e-03 want= 1.01000000e-03 OK axial: reaction0x got=-1.00000000e+03 want=-1.00000000e+03 OK axial: util got= 4.00000000e-03 want= 4.00000000e-03 OK prestrain: N (compression<0) got=-6.66666667e+03 want=-6.66666667e+03 OK prestrain: u1x (+x) got= 6.66666667e-05 want= 6.66666667e-05 OK prestrain: u0x (-x) got=-6.66666667e-05 want=-6.66666667e-05 OK thermal: N got=-6.66666667e+03 want=-6.66666667e+03 OK settle: u1x got= 5.00000250e-04 want= 5.00000000e-04 OK settle: N got= 4.99999750e+04 want= 5.00000000e+04 OK settle: reaction1x got= 4.99999750e+04 want= 5.00000000e+04 OK cable slack: active got= 0.00000000e+00 want= 0.00000000e+00 OK cable slack: force=0 got= 0.00000000e+00 want= 0.00000000e+00 OK cable slack: util=0 got= 0.00000000e+00 want= 0.00000000e+00 OK cable slack: u1x got=-1.00000000e-03 want=-1.00000000e-03 OK cable taut: active got= 1.00000000e+00 want= 1.00000000e+00 OK cable taut: N got= 1.00000000e+03 want= 1.00000000e+03 OK geom: N (tension) got= 1.99999600e+05 want= 2.00000000e+05 OK geom: u1y w/ N/L stiffening got= 9.09092563e-03 want= 9.09092562e-03 OK ALL PASS (0 failures) [stdout] axial: N got= 1.00000000e+03 want= 1.00000000e+03 OK axial: u0x got= 1.00000000e-03 want= 1.00000000e-03 OK axial: u1x got= 1.01000000e-03 want= 1.01000000e-03 OK axial: reaction0x got=-1.00000000e+03 want=-1.00000000e+03 OK axial: util got= 4.00000000e-03 want= 4.00000000e-03 OK prestrain: N (compression<0) got=-6.66666667e+03 want=-6.66666667e+03 OK prestrain: u1x (+x) got= 6.66666667e-05 want= 6.66666667e-05 OK prestrain: u0x (-x) got=-6.66666667e-05 want=-6.66666667e-05 OK thermal: N got=-6.66666667e+03 want=-6.66666667e+03 OK settle: u1x got= 5.00000250e-04 want= 5.00000000e-04 OK settle: N got= 4.99999750e+04 want= 5.00000000e+04 OK settle: reaction1x got= 4.99999750e+04 want= 5.00000000e+04 OK cable slack: active got= 0.00000000e+00 want= 0.00000000e+00 OK cable slack: force=0 got= 0.00000000e+00 want= 0.00000000e+00 OK cable slack: util=0 got= 0.00000000e+00 want= 0.00000000e+00 OK cable slack: u1x got=-1.00000000e-03 want=-1.00000000e-03 OK cable taut: active got= 1.00000000e+00 want= 1.00000000e+00 OK cable taut: N got= 1.00000000e+03 want= 1.00000000e+03 OK geom: N (tension) got= 1.99999600e+05 want= 2.00000000e+05 OK geom: u1y w/ N/L stiffening got= 9.09092563e-03 want= 9.09092562e-03 OK ALL PASS (0 failures)
cat > /tmp/verify2.cpp <<'EOF'
#include <cmath>
#include <cstdio>
#include <string>
#include <vector>
#include "truss2d/model.hpp"
#include "truss2d/solver.hpp"
#include "truss2d/dense_matrix.hpp"
using namespace truss2d;
static int failures=0;
static void chk(const std::string&n,double g,double w,double tol){
double e=std::fabs(g-w); double rel=e/(std::fabs(w)+1e-30);
bool ok=e<=tol||rel<=1e-6; if(!ok)++failures;
std::printf("%-42s got=% .6e want=% .6e %s\n",n.c_str(),g,w,ok?"OK":" <<<FAIL");}
int main(){
const double E=2e11,A=1e-3;
// ---- 2D symmetric two-bar truss, statics ----
{
double P=-1000;
Model m;
m.add_node(-1,1); m.add_node(1,1); m.add_node(0,0);
m.add_element(0,2,Section{E,A});
m.add_element(1,2,Section{E,A});
m.add_spring(0,1,0,1e14); m.add_spring(0,0,1,1e14);
m.add_spring(1,1,0,1e14); m.add_spring(1,0,1,1e14);
m.add_load(5,P); // node2 y
auto r=StaticSolver(m).solve();
double Nexp=1000.0/std::sqrt(2.0);
chk("2Dtruss: N0",r.axial_forces[0],Nexp,1e-2);
chk("2Dtruss: N1",r.axial_forces[1],Nexp,1e-2);
chk("2Dtruss: node2 x=0",r.displacements[4],0.0,1e-9);
}
// ---- skew roller (inclined spring) ----
{
// node0 free, single inclined spring at 45 deg, load along x.
// spring resists only along dir=(1,1)/sqrt2. Load Px at node0.
// Equilibrium: kappa*(u.dir)dir = F. Only the dir-component is resisted;
// perpendicular is unconstrained -> singular unless load is along dir.
// Use two nonparallel springs => 2D support.
double kx=2e6, P=1000;
Model m;
m.add_node(0,0); m.add_node(1,0);
m.add_element(0,1,Section{E,A});
m.add_spring(0,1,1,kx); // 45 deg skew
m.add_spring(0,1,-1,kx); // -45 deg -> together act 2D
m.add_spring(1,0,1,1e14);
m.add_load(2,P); // node1 x
auto r=StaticSolver(m).solve();
// node0 held by two springs each kappa along (1,±1)/sqrt2:
// sum kappa dir dir^T = kx*( [.5 .5;.5 .5]+[.5 -.5;-.5 .5] )=kx*[1 0;0 1].
// so effective 2D stiffness kx in x,y. Bar to node1.
// x-system: node0x (kx + EA/L) - EA/L node1x =0 ; node1x: -EA/L n0 + EA/L n1 = P
double EA=E*A,L=1.0;
// solve 2x2
double a=kx+EA/L,b=-EA/L,c=-EA/L,d=EA/L,det=a*d-b*c;
double u0=( d*0 - b*P)/det, u1=(-c*0 + a*P)/det;
chk("skew: u0x",r.displacements[0],u0,1e-9);
chk("skew: u1x",r.displacements[2],u1,1e-9);
chk("skew: N",r.axial_forces[0],EA*(u1-u0)/L,1e-3);
chk("skew: node0 y=0",r.displacements[1],0.0,1e-9);
}
// ---- strut: stays in compression, drops in tension ----
{
double L=2.0;
Model m;
m.add_node(0,0); m.add_node(L,0);
m.add_element(0,1,Section{E,A},MemberKind::Strut);
m.add_spring(0,1,0,1e14); m.add_spring(0,0,1,1e14);
m.add_spring(1,1,0,1e6); m.add_spring(1,0,1,1e14);
m.add_load(2,-1000); // push node1 toward node0 -> compression
auto r=StaticSolver(m).solve();
chk("strut compress: active",r.active[0]?1:0,1,0);
chk("strut compress: N<0",r.axial_forces[0],-1000.0,1e-2);
}
{
double L=2.0;
Model m;
m.add_node(0,0); m.add_node(L,0);
m.add_element(0,1,Section{E,A},MemberKind::Strut);
m.add_spring(0,1,0,1e14); m.add_spring(0,0,1,1e14);
m.add_spring(1,1,0,1e6); m.add_spring(1,0,1,1e14);
m.add_load(2,1000); // pull node1 away -> would be tension -> strut separates
auto r=StaticSolver(m).solve();
chk("strut tension: inactive",r.active[0]?1:0,0,0);
chk("strut tension: N=0",r.axial_forces[0],0.0,0);
chk("strut tension: u1x",r.displacements[2],1000.0/1e6,1e-12);
}
// ---- equilibrium residual via public helpers (no geometric, no unilat) ----
{
Model m;
m.add_node(0,0); m.add_node(3,0); m.add_node(3,4);
m.add_element(0,1,Section{E,A});
m.add_element(1,2,Section{E,A});
m.add_element(0,2,Section{E,A});
m.add_spring(0,1,0,1e13); m.add_spring(0,0,1,1e13);
m.add_spring(1,0,1,1e13); m.add_spring(2,1,0,1e13);
m.add_load(5,-2000);
StaticSolver s(m);
auto r=s.solve();
std::vector<bool> act=r.active;
DenseMatrix K=s.assemble_system(act);
std::vector<double> fp=s.prestrain_force(act);
// residual = K u - (loads + prestrain + settle) - (reactions)
// Since elastic-only (no prestrain, no geometric), K u should equal
// loads + reactions. Check K u - loads - reactions ~ 0.
std::size_t n=m.num_dofs();
std::vector<double> Ku(n,0);
for(std::size_t i=0;i<n;i++){double sgn=0;for(std::size_t j=0;j<n;j++)sgn+=K(i,j)*r.displacements[j];Ku[i]=sgn;}
std::vector<double> load(n,0); for(auto&ld:m.loads())load[ld.dof]+=ld.value;
double maxres=0;
for(std::size_t i=0;i<n;i++){double res=Ku[i]-fp[i]-load[i]-r.reactions[i];maxres=std::max(maxres,std::fabs(res));}
chk("residual: K u - f - R ~ 0",maxres,0.0,1e-3);
}
std::printf("\n%s (%d failures)\n",failures?"SOME FAIL":"ALL PASS",failures);
return failures?1:0;
}
EOF
g++ -std=c++17 -Iinclude /tmp/verify2.cpp src/model.cpp src/dense_matrix.cpp src/solver.cpp -o /tmp/verify2 2>&1 | head && /tmp/verify2Exit code 1 2Dtruss: N0 got= 7.071043e+02 want= 7.071068e+02 OK 2Dtruss: N1 got= 7.071043e+02 want= 7.071068e+02 OK 2Dtruss: node2 x=0 got= 0.000000e+00 want= 0.000000e+00 OK skew: u0x got= 5.000000e-04 want= 5.000000e-04 OK skew: u1x got= 5.050000e-04 want= 5.050000e-04 OK skew: N got= 1.000000e+03 want= 1.000000e+03 OK skew: node0 y=0 got= 0.000000e+00 want= 0.000000e+00 OK strut compress: active got= 1.000000e+00 want= 1.000000e+00 OK strut compress: N<0 got=-9.900990e+02 want=-1.000000e+03 <<<FAIL strut tension: inactive got= 0.000000e+00 want= 0.000000e+00 OK strut tension: N=0 got= 0.000000e+00 want= 0.000000e+00 OK strut tension: u1x got= 1.000000e-03 want= 1.000000e-03 OK residual: K u - f - R ~ 0 got= 1.322753e+03 want= 0.000000e+00 <<<FAIL SOME FAIL (2 failures) [error] tool reported failure
cat > /tmp/verify3.cpp <<'EOF'
#include <cmath>
#include <cstdio>
#include <vector>
#include "truss2d/model.hpp"
#include "truss2d/solver.hpp"
#include "truss2d/dense_matrix.hpp"
using namespace truss2d;
int main(){
const double E=2e11,A=1e-3;
Model m;
m.add_node(0,0); m.add_node(3,0); m.add_node(3,4); m.add_node(6,0);
// mix: a prestrained bar, a normal bar, a cable
m.add_element(0,1,Section{E,A},MemberKind::Bar,-3e-4); // pretension
m.add_element(1,2,Section{E,A},MemberKind::Bar);
m.add_element(0,2,Section{E,A},MemberKind::Cable);
m.add_element(1,3,Section{E,A},MemberKind::Bar,0,1.2e-5,-30); // thermal
m.add_spring(0,1,0,1e13); m.add_spring(0,0,1,1e13);
m.add_spring(2,0,1,1e13); m.add_spring(3,1,0,1e13); m.add_spring(3,0,1,1e13);
m.add_spring(2,1,0,5e6); // skew finite
m.add_load(2,-1500); m.add_load(5,-2500);
StaticSolver s(m);
auto r=s.solve();
std::size_t n=m.num_dofs(), ne=m.num_elements();
std::vector<bool> act=r.active;
// elastic-only stiffness K_e = sum active element_stiffness_global
DenseMatrix Ke(n,n);
for(std::size_t e=0;e<ne;e++){ if(!act[e])continue;
DenseMatrix k=s.element_stiffness_global(e); const auto&el=m.elements()[e];
std::size_t g[4]={2*el.n1,2*el.n1+1,2*el.n2,2*el.n2+1};
for(int i=0;i<4;i++)for(int j=0;j<4;j++)Ke(g[i],g[j])+=k(i,j);}
std::vector<double> fp=s.prestrain_force(act);
std::vector<double> load(n,0); for(auto&ld:m.loads())load[ld.dof]+=ld.value;
// nodal equilibrium: K_e u - fp - load - R = 0
double maxres=0; std::size_t worst=0;
for(std::size_t i=0;i<n;i++){
double Ku=0; for(std::size_t j=0;j<n;j++)Ku+=Ke(i,j)*r.displacements[j];
double res=Ku-fp[i]-load[i]-r.reactions[i];
if(std::fabs(res)>maxres){maxres=std::fabs(res);worst=i;}
}
std::printf("max nodal-equilibrium residual = %.3e at dof %zu\n",maxres,worst);
// also verify axial forces recovered match N=EA(delta/L - eps0) and util sign
std::printf("forces: ");
for(std::size_t e=0;e<ne;e++) std::printf("[%zu]N=%.3e act=%d util=%.3e ",e,r.axial_forces[e],(int)r.active[e],r.utilization[e]);
std::printf("\n%s\n", maxres<1e-2?"EQUILIBRIUM OK":"EQUILIBRIUM FAIL");
return maxres<1e-2?0:1;
}
EOF
g++ -std=c++17 -Iinclude /tmp/verify3.cpp src/model.cpp src/dense_matrix.cpp src/solver.cpp -o /tmp/verify3 2>&1 | head && /tmp/verify3max nodal-equilibrium residual = 1.099e-05 at dof 3 forces: [0]N=6.525e+04 act=1 util=2.610e-01 [1]N=-1.099e-05 act=1 util=-4.396e-11 [2]N=0.000e+00 act=0 util=0.000e+00 [3]N=6.675e+04 act=1 util=2.670e-01 EQUILIBRIUM OK [stdout] max nodal-equilibrium residual = 1.099e-05 at dof 3 forces: [0]N=6.525e+04 act=1 util=2.610e-01 [1]N=-1.099e-05 act=1 util=-4.396e-11 [2]N=0.000e+00 act=0 util=0.000e+00 [3]N=6.675e+04 act=1 util=2.670e-01 EQUILIBRIUM OK
rm -f verify.cpp && cmake --build build --parallel >/tmp/b.log 2>&1 && echo "BUILD OK" && grep -E "warning|error" /tmp/b.log | head; cd build && ctest --output-on-failure 2>&1 | tail -5
BUILD OK 1/1 Test #1: smoke ............................ Passed 0.00 sec 100% tests passed, 0 tests failed out of 1 Total Test time (real) = 0.00 sec [stdout] BUILD OK 1/1 Test #1: smoke ............................ Passed 0.00 sec 100% tests passed, 0 tests failed out of 1 Total Test time (real) = 0.00 sec
No file changes, the agent left the workspace untouched.
== Restoring pristine project ==
== Staging candidate solver ==
== Injecting hidden grading tests ==
== Configuring (cmake) ==
-- The CXX compiler identification is GNU 11.4.0
-- Detecting CXX compiler ABI info
-- Detecting CXX compiler ABI info - done
-- Check for working CXX compiler: /usr/bin/c++ - skipped
-- Detecting CXX compile features
-- Detecting CXX compile features - done
-- Configuring done
-- Generating done
-- Build files have been written to: /tmp/tmp.aMsqVN2jTq/truss2d/build_grade
== Building ==
[ 5%] Building CXX object CMakeFiles/truss2d.dir/src/dense_matrix.cpp.o
[ 10%] Building CXX object CMakeFiles/truss2d.dir/src/solver.cpp.o
[ 15%] Building CXX object CMakeFiles/truss2d.dir/src/io.cpp.o
[ 21%] Building CXX object CMakeFiles/truss2d.dir/src/model.cpp.o
[ 26%] Linking CXX static library libtruss2d.a
[ 26%] Built target truss2d
[ 31%] Building CXX object tests/CMakeFiles/test_soak.dir/test_soak.cpp.o
[ 36%] Building CXX object tests/CMakeFiles/test_unilateral.dir/test_unilateral.cpp.o
[ 42%] Building CXX object tests/CMakeFiles/test_geometric.dir/test_geometric.cpp.o
[ 47%] Building CXX object tests/CMakeFiles/test_core.dir/test_core.cpp.o
[ 52%] Building CXX object CMakeFiles/truss2d_cli.dir/src/main.cpp.o
[ 57%] Building CXX object tests/CMakeFiles/test_degenerate.dir/test_degenerate.cpp.o
[ 63%] Building CXX object tests/CMakeFiles/test_settlement.dir/test_settlement.cpp.o
[ 68%] Linking CXX executable truss2d_cli
[ 68%] Built target truss2d_cli
[ 73%] Linking CXX executable test_degenerate
[ 73%] Built target test_degenerate
[ 78%] Linking CXX executable test_settlement
[ 84%] Linking CXX executable test_soak
[ 89%] Linking CXX executable test_core
[ 89%] Built target test_soak
[ 89%] Built target test_settlement
[ 89%] Built target test_core
[ 94%] Linking CXX executable test_unilateral
[ 94%] Built target test_unilateral
[100%] Linking CXX executable test_geometric
[100%] Built target test_geometric
== Running hidden tests ==
Test project /tmp/tmp.aMsqVN2jTq/truss2d/build_grade
Start 1: test_core
1/6 Test #1: test_core ........................ Passed 0.00 sec
Start 2: test_unilateral
2/6 Test #2: test_unilateral .................. Passed 0.00 sec
Start 3: test_geometric
3/6 Test #3: test_geometric ...................***Failed 0.00 sec
[ FAIL ] geometric_pretensioned_lateral_stiffness_vs_analytic: unexpected exception: solve_spd: matrix is singular or not positive-definite
[ FAIL ] geometric_uses_current_total_force: unexpected exception: solve_spd: matrix is singular or not positive-definite
[ PASS ] geometric_assemble_system_remains_elastic_only
[ FAIL ] geometric_force_recovery_is_axial_only: unexpected exception: solve_spd: matrix is singular or not positive-definite
[ FAIL ] geometric_offaxis_pretensioned_net_vs_analytic: unexpected exception: solve_spd: matrix is singular or not positive-definite
[ FAIL ] geometric_multi_segment_chain_axial_vs_oracle: unexpected exception: solve_spd: matrix is singular or not positive-definite
----
1/6 tests passed
Start 4: test_soak
4/6 Test #4: test_soak ........................***Failed 0.00 sec
[ FAIL ] soak_pretensioned_thermal_cable_nets: unexpected exception: solve_spd: matrix is singular or not positive-definite
----
0/1 tests passed
Start 5: test_degenerate
5/6 Test #5: test_degenerate .................. Passed 0.00 sec
Start 6: test_settlement
6/6 Test #6: test_settlement .................. Passed 0.00 sec
67% tests passed, 2 tests failed out of 6
Label Time Summary:
hidden = 0.01 sec*proc (6 tests)
Total Test time (real) = 0.02 sec
The following tests FAILED:
3 - test_geometric (Failed)
4 - test_soak (Failed)
Errors while running CTest
FAIL: hidden tests failedReproduce this trial: git checkout 2f94510 && PYTHONPATH=src python3 scripts/build_site.py , then open trial/trial_82dc1bf4532a4fda. Re-running the agent live requires EVAL_PLATFORM_ENABLE_OAUTH_SMOKE=1 and is non-deterministic.
Trial trial_82dc1bf4532a4fda · verifier authoritative; classifier explanatory.