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/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/instruction.md
File does not exist. Note: your current working directory is /opt/truss2d. [error] tool reported failure
**/*test*
tests/test_util.hpp tests/test_smoke.cpp
/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
ls -R /opt/truss2d --ignore=build 2>/dev/null; echo "---IO---"; cat /opt/truss2d/src/io.cpp 2>/dev/null | head -200
/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
---IO---
#include "truss2d/io.hpp"
#include <fstream>
#include <iomanip>
#include <iostream>
#include <sstream>
#include <stdexcept>
#include <string>
namespace truss2d {
namespace {
[[noreturn]] void fail(std::size_t line_no, const std::string& msg) {
std::ostringstream oss;
oss << "parse_model: line " << line_no << ": " << msg;
throw std::runtime_error(oss.str());
}
const char* kind_name(MemberKind k) {
switch (k) {
case MemberKind::Bar: return "bar";
case MemberKind::Cable: return "cable";
case MemberKind::Strut: return "strut";
}
return "?";
}
} // namespace
Model parse_model(std::istream& in) {
Model model;
std::string line;
std::size_t line_no = 0;
while (std::getline(in, line)) {
++line_no;
const auto hash = line.find('#');
if (hash != std::string::npos) {
line.erase(hash);
}
std::istringstream ls(line);
std::string tag;
if (!(ls >> tag)) {
continue; // blank / comment-only line
}
if (tag == "node") {
double x, y;
if (!(ls >> x >> y)) fail(line_no, "expected: node <x> <y>");
model.add_node(x, y);
} else if (tag == "elem") {
std::size_t n1, n2;
double E, A;
if (!(ls >> n1 >> n2 >> E >> A)) {
fail(line_no,
"expected: elem <n1> <n2> <E> <A> [kind] [prestrain] [alpha] [dT]");
}
int kind_code = 0;
double prestrain = 0.0;
double alpha = 0.0;
double dT = 0.0;
ls >> kind_code; // optional
ls >> prestrain; // optional
ls >> alpha; // optional
ls >> dT; // optional
MemberKind kind = MemberKind::Bar;
if (kind_code == 1) kind = MemberKind::Cable;
else if (kind_code == 2) kind = MemberKind::Strut;
else if (kind_code != 0) fail(line_no, "kind must be 0/1/2 (bar/cable/strut)");
try {
model.add_element(n1, n2, Section{E, A}, kind, prestrain, alpha, dT);
} catch (const std::exception& ex) {
fail(line_no, ex.what());
}
} else if (tag == "spring") {
std::size_t node;
double dx, dy, kappa;
if (!(ls >> node >> dx >> dy >> kappa)) {
fail(line_no, "expected: spring <node> <dx> <dy> <kappa>");
}
try {
model.add_spring(node, dx, dy, kappa);
} catch (const std::exception& ex) {
fail(line_no, ex.what());
}
} else if (tag == "load") {
std::size_t node, comp;
double value;
if (!(ls >> node >> comp >> value)) {
fail(line_no, "expected: load <node> <comp> <value>");
}
if (comp > 1) fail(line_no, "component must be 0 (x) or 1 (y)");
try {
model.add_load(2 * node + comp, value);
} catch (const std::exception& ex) {
fail(line_no, ex.what());
}
} else {
fail(line_no, "unknown record '" + tag + "'");
}
}
return model;
}
Model parse_model_file(const std::string& path) {
std::ifstream in(path);
if (!in) {
throw std::runtime_error("parse_model_file: cannot open '" + path + "'");
}
return parse_model(in);
}
void write_report(std::ostream& out, const Model& model, const SolveResult& result) {
out << std::scientific << std::setprecision(6);
out << "# Truss2D static analysis report\n";
out << "nodes=" << model.num_nodes()
<< " elements=" << model.num_elements()
<< " dofs=" << model.num_dofs() << "\n\n";
out << "[displacements] (node, ux, uy)\n";
for (std::size_t i = 0; i < model.num_nodes(); ++i) {
out << i << " " << result.displacements[2 * i] << " "
<< result.displacements[2 * i + 1] << "\n";
}
out << "\n[reactions] (node, Rx, Ry)\n";
for (std::size_t i = 0; i < model.num_nodes(); ++i) {
out << i << " " << result.reactions[2 * i] << " "
<< result.reactions[2 * i + 1] << "\n";
}
out << "\n[axial_forces] (element, kind, active, N [+tension], utilization)\n";
for (std::size_t e = 0; e < model.num_elements(); ++e) {
const bool act = e < result.active.size() ? result.active[e] : true;
const double util = e < result.utilization.size() ? result.utilization[e] : 0.0;
out << e << " " << kind_name(model.elements()[e].kind) << " "
<< (act ? 1 : 0) << " " << result.axial_forces[e] << " "
<< util << "\n";
}
}
} // namespace truss2d
[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
---IO---
#include "truss2d/io.hpp"
#include <fstream>
#include <iomanip>
#include <iostream>
#include <sstream>
#include <stdexcept>
#include <string>
namespace truss2d {
namespace {
[[noreturn]] void fail(std::size_t line_no, const std::string& msg) {
std::ostringstream oss;
oss << "parse_model: line " << line_no << ": " << msg;
throw std::runtime_error(oss.str());
}
const char* kind_name(MemberKind k) {
switch (k) {
case MemberKind::Bar: return "bar";
case MemberKind::Cable: return "cable";
case MemberKind::Strut: return "strut";
}
return "?";
}
} // namespace
Model parse_model(std::istream& in) {
Model model;
std::string line;
std::size_t line_no = 0;
while (std::getline(in, line)) {
++line_no;
const auto hash = line.find('#');
if (hash != std::string::npos) {
line.erase(hash);
}
std::istringstream ls(line);
std::string tag;
if (!(ls >> tag)) {
continue; // blank / comment-only line
}
if (tag == "node") {
double x, y;
if (!(ls >> x >> y)) fail(line_no, "expected: node <x> <y>");
model.add_node(x, y);
} else if (tag == "elem") {
std::size_t n1, n2;
double E, A;
if (!(ls >> n1 >> n2 >> E >> A)) {
fail(line_no,
"expected: elem <n1> <n2> <E> <A> [kind] [prestrain] [alpha] [dT]");
}
int kind_code = 0;
double prestrain = 0.0;
double alpha = 0.0;
double dT = 0.0;
ls >> kind_code; // optional
ls >> prestrain; // optional
ls >> alpha; // optional
ls >> dT; // optional
MemberKind kind = MemberKind::Bar;
if (kind_code == 1) kind = MemberKind::Cable;
else if (kind_code == 2) kind = MemberKind::Strut;
else if (kind_code != 0) fail(line_no, "kind must be 0/1/2 (bar/cable/strut)");
try {
model.add_element(n1, n2, Section{E, A}, kind, prestrain, alpha, dT);
} catch (const std::exception& ex) {
fail(line_no, ex.what());
}
} else if (tag == "spring") {
std::size_t node;
double dx, dy, kappa;
if (!(ls >> node >> dx >> dy >> kappa)) {
fail(line_no, "expected: spring <node> <dx> <dy> <kappa>");
}
try {
model.add_spring(node, dx, dy, kappa);
} catch (const std::exception& ex) {
fail(line_no, ex.what());
}
} else if (tag == "load") {
std::size_t node, comp;
double value;
if (!(ls >> node >> comp >> value)) {
fail(line_no, "expected: load <node> <comp> <value>");
}
if (comp > 1) fail(line_no, "component must be 0 (x) or 1 (y)");
try {
model.add_load(2 * node + comp, value);
} catch (const std::exception& ex) {
fail(line_no, ex.what());
}
} else {
fail(line_no, "unknown record '" + tag + "'");
}
}
return model;
}
Model parse_model_file(const std::string& path) {
std::ifstream in(path);
if (!in) {
throw std::runtime_error("parse_model_file: cannot open '" + path + "'");
}
return parse_model(in);
}
void write_report(std::ostream& out, const Model& model, const SolveResult& result) {
out << std::scientific << std::setprecision(6);
out << "# Truss2D static analysis report\n";
out << "nodes=" << model.num_nodes()
<< " elements=" << model.num_elements()
<< " dofs=" << model.num_dofs() << "\n\n";
out << "[displacements] (node, ux, uy)\n";
for (std::size_t i = 0; i < model.num_nodes(); ++i) {
out << i << " " << result.displacements[2 * i] << " "
<< result.displacements[2 * i + 1] << "\n";
}
out << "\n[reactions] (node, Rx, Ry)\n";
for (std::size_t i = 0; i < model.num_nodes(); ++i) {
out << i << " " << result.reactions[2 * i] << " "
<< result.reactions[2 * i + 1] << "\n";
}
out << "\n[axial_forces] (element, kind, active, N [+tension], utilization)\n";
for (std::size_t e = 0; e < model.num_elements(); ++e) {
const bool act = e < result.active.size() ? result.active[e] : true;
const double util = e < result.utilization.size() ? result.utilization[e] : 0.0;
out << e << " " << kind_name(model.elements()[e].kind) << " "
<< (act ? 1 : 0) << " " << result.axial_forces[e] << " "
<< util << "\n";
}
}
} // namespace truss2d---
name: truss2d-solver-model
description: Physics/sign conventions for the truss2d StaticSolver implementation
metadata:
type: project
---
truss2d static solver (`src/solver.cpp`) implements a 2D two-force-member network with cables/struts/bars, spring supports, prestrain, and initial-stress (geometric) stiffness.
Key conventions derived from the ticket:
- Local DOF order {n1x,n1y,n2x,n2y}. Direction (c,s) from n1→n2. b=(-c,-s,c,s), δ=bᵀu.
- Elastic 4x4: k·r·rᵀ, r=(c,s,-c,-s), k=EA/L.
- Installed stress-free strain ε_free = prestrain + alpha·dT.
- Prestrain force (RHS): EA·ε_free·b. Recovered N = k·δ − EA·ε_free (tension +).
Too-long/heated (ε_free>0) → compression; too-short → tension.
- Geometric (initial-stress) stiffness: (N/L)·a·aᵀ, a=(s,-c,-s,c). Only active members, uses total N. ONLY in solve(), NOT in helpers.
- Spring: stiffness κ·d·dᵀ at node DOFs. Settlement RHS force +κ·s0·d. Reaction = −κ(dᵀu_node − s0)·d.
- solve() iterates active-set (cable active iff trial N>0, strut iff N<0, bar always) + fixed-point on geometric N until both stable.
- Utilization = N/(2.5e8·A), signed. Inactive unilateral → 0 force, 0 util.
- solve() must call model_.validate() first; let solve_spd throw on non-SPD.
test -f /opt/truss2d/MEMORY.md && echo exists || echo none; cat /logs/agent/sessions/projects/-opt-truss2d/memory/MEMORY.md 2>/dev/null || echo "no memory index"
none no memory index [stdout] none no memory index
- [Truss2D solver model](truss2d-solver-model.md) , physics & sign conventions for src/solver.cpp
#include "truss2d/solver.hpp"
#include <algorithm>
#include <cmath>
#include <cstddef>
#include <vector>
#include "truss2d/dense_matrix.hpp"
#include "truss2d/types.hpp"
namespace truss2d {
namespace {
// Per-element geometry/material cache used by solve().
//
// Local DOF order is {n1.x, n1.y, n2.x, n2.y}. The axial extension produced by
// nodal displacement is delta = b . u_local with b = (-c, -s, c, s); the
// ordinary elastic 4x4 is k * b b^T (equivalently k * r r^T with r = -b).
struct ElemGeom {
std::size_t d[4]; // global DOF indices of the 4 local DOFs
double c; // direction cosine (n1 -> n2)
double s; // direction sine
double L; // undeformed length
double EA; // E * A
double k; // EA / L
double eps_free; // installed stress-free strain (mechanical + thermal)
};
} // namespace
DenseMatrix StaticSolver::element_stiffness_global(std::size_t e) const {
const Vec2 dir = model_.element_direction(e);
const Element& el = model_.elements()[e];
const double k = el.section.E * el.section.A / model_.element_length(e);
// K_e = k * r r^T with r = (c, s, -c, -s).
const double r[4] = {dir.x, dir.y, -dir.x, -dir.y};
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 * r[i] * r[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);
// Ordinary elastic stiffness of every participating member.
const auto& elems = model_.elements();
for (std::size_t e = 0; e < elems.size(); ++e) {
if (e < active.size() && !active[e]) continue;
const DenseMatrix Ke = element_stiffness_global(e);
const std::size_t g[4] = {2 * elems[e].n1, 2 * elems[e].n1 + 1,
2 * elems[e].n2, 2 * elems[e].n2 + 1};
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);
}
}
}
// Grounded spring supports: kappa * d d^T on the node's two DOFs.
for (const auto& sp : model_.springs()) {
const std::size_t a = 2 * sp.node;
const std::size_t b = 2 * sp.node + 1;
const double dx = sp.dir.x;
const double dy = sp.dir.y;
K(a, a) += sp.kappa * dx * dx;
K(a, b) += sp.kappa * dx * dy;
K(b, a) += sp.kappa * dy * dx;
K(b, b) += sp.kappa * dy * dy;
}
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) {
if (e < active.size() && !active[e]) continue;
const Element& el = elems[e];
const Vec2 dir = model_.element_direction(e);
const double EA = el.section.E * el.section.A;
const double eps_free = el.prestrain + el.alpha * el.dT;
const double mag = EA * eps_free; // F = EA*eps_free * b
const std::size_t g[4] = {2 * el.n1, 2 * el.n1 + 1, 2 * el.n2,
2 * el.n2 + 1};
const double b[4] = {-dir.x, -dir.y, dir.x, dir.y};
for (std::size_t i = 0; i < 4; ++i) {
f[g[i]] += mag * b[i];
}
}
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();
// Cache per-element geometry/material.
std::vector<ElemGeom> eg(ne);
for (std::size_t e = 0; e < ne; ++e) {
const Element& el = elems[e];
const Vec2 dir = model_.element_direction(e);
ElemGeom g;
g.d[0] = 2 * el.n1;
g.d[1] = 2 * el.n1 + 1;
g.d[2] = 2 * el.n2;
g.d[3] = 2 * el.n2 + 1;
g.c = dir.x;
g.s = dir.y;
g.L = model_.element_length(e);
g.EA = el.section.E * el.section.A;
g.k = g.EA / g.L;
g.eps_free = el.prestrain + el.alpha * el.dT;
eg[e] = g;
}
// External nodal load vector.
std::vector<double> Fext(n, 0.0);
for (const auto& ld : model_.loads()) {
Fext[ld.dof] += ld.value;
}
// Constant settlement contribution to the right-hand side:
// +kappa * settlement * dir at the supported node's DOFs.
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;
}
// Recover the total axial force of element e (tension positive) from a
// displacement field: N = k * delta - EA * eps_free.
auto recover = [&](std::size_t e, const std::vector<double>& u) {
const ElemGeom& g = eg[e];
const double delta = -g.c * u[g.d[0]] - g.s * u[g.d[1]] +
g.c * u[g.d[2]] + g.s * u[g.d[3]];
return g.k * delta - g.EA * g.eps_free;
};
std::vector<bool> active(ne, true);
std::vector<double> N(ne, 0.0); // member force used for the geometric term
std::vector<double> u(n, 0.0);
// Resolve the active set and the force-dependent initial-stress effect
// together: the assumed forces, the displacements they produce, and the
// forces recovered from those displacements must all coincide.
const int kMaxIter = 1000;
for (int it = 0; it < kMaxIter; ++it) {
// Ordinary elastic + spring stiffness for the current active set ...
DenseMatrix K = assemble_system(active);
// ... plus the initial-stress (geometric) stiffness of active members,
// (N/L) * a a^T with a = (s, -c, -s, c) acting on transverse motion.
for (std::size_t e = 0; e < ne; ++e) {
if (!active[e]) continue;
const ElemGeom& g = eg[e];
const double coef = N[e] / g.L;
if (coef == 0.0) continue;
const double a[4] = {g.s, -g.c, -g.s, g.c};
for (std::size_t i = 0; i < 4; ++i) {
for (std::size_t j = 0; j < 4; ++j) {
K(g.d[i], g.d[j]) += coef * a[i] * a[j];
}
}
}
std::vector<double> F(n);
const std::vector<double> fp = prestrain_force(active);
for (std::size_t i = 0; i < n; ++i) {
F[i] = Fext[i] + fp[i] + Fsettle[i];
}
u = solve_spd(K, F);
// Trial forces from this field, and the implied active set.
std::vector<bool> newactive(ne, true);
std::vector<double> Nnext(ne, 0.0);
double maxdiff = 0.0;
double scale = 1.0;
for (std::size_t e = 0; e < ne; ++e) {
const double Ntrial = recover(e, u);
bool act = true;
switch (elems[e].kind) {
case MemberKind::Bar: act = true; break;
case MemberKind::Cable: act = (Ntrial > 0.0); break;
case MemberKind::Strut: act = (Ntrial < 0.0); break;
}
newactive[e] = act;
Nnext[e] = act ? Ntrial : 0.0;
maxdiff = std::max(maxdiff, std::fabs(Nnext[e] - N[e]));
scale = std::max(scale, std::fabs(Nnext[e]));
}
const bool active_stable = (newactive == active);
active = newactive;
N = Nnext;
if (active_stable && maxdiff <= 1e-12 * scale) break;
}
// Assemble outputs.
SolveResult result;
result.displacements = u;
result.reactions.assign(n, 0.0);
for (const auto& sp : model_.springs()) {
const double along =
sp.dir.x * u[2 * sp.node] + sp.dir.y * u[2 * sp.node + 1];
const double eta = along - sp.settlement; // stored stretch
const double R = -sp.kappa * eta; // force on structure along dir
result.reactions[2 * sp.node] += R * sp.dir.x;
result.reactions[2 * sp.node + 1] += R * sp.dir.y;
}
result.axial_forces.assign(ne, 0.0);
result.utilization.assign(ne, 0.0);
result.active.assign(ne, true);
const double sigma_allow = 2.5e8;
for (std::size_t e = 0; e < ne; ++e) {
result.active[e] = active[e];
if (active[e]) {
result.axial_forces[e] = N[e];
const double capacity = sigma_allow * elems[e].section.A;
result.utilization[e] = N[e] / capacity;
}
}
return result;
}
} // namespace truss2d
cmake -S . -B build -DCMAKE_BUILD_TYPE=Release >/tmp/cm.log 2>&1 && cmake --build build --parallel >/tmp/build.log 2>&1 && echo BUILD_OK && cd build && ctest --output-on-failure 2>&1 | tail -30 || { echo FAILED; tail -40 /tmp/build.log; }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 seccat > /tmp/verify.cpp <<'EOF'
#include <cmath>
#include <cstdio>
#include <vector>
#include "truss2d/model.hpp"
#include "truss2d/solver.hpp"
using namespace truss2d;
int fails=0;
void chk(bool c,const char*m){ if(!c){printf("FAIL: %s\n",m);++fails;} else printf("ok: %s\n",m);}
void near(double g,double w,double tol,const char*m){ if(std::fabs(g-w)>tol){printf("FAIL: %s got %.10g want %.10g\n",m,g,w);++fails;} else printf("ok: %s (%.6g)\n",m,g);}
// residual = K(active)*u - prestrain(active) - Fsettle - Fext ; should be ~ -reactions?
int main(){
double E=2e11, A=1e-3, EA=E*A;
// --- Case 1: single horizontal bar, node0 fully sprung, node1 y-sprung, load +x at node1 ---
{
Model m; m.add_node(0,0); m.add_node(2,0);
m.add_element(0,1,Section{E,A});
double K=1e15; // very stiff springs ~ fixed
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 L=2.0, k=EA/L;
// node1 x disp ~ P/k (springs near rigid). axial = k*u = P (tension)
near(r.axial_forces[0], P, 1.0, "C1 axial=P tension");
near(r.displacements[2], P/k, 1e-6, "C1 u1x = P/k");
near(r.reactions[0], -P, 1.0, "C1 Rx node0 = -P");
chk(r.active[0],"C1 bar active");
near(r.utilization[0], P/(2.5e8*A), 1e-9,"C1 util");
}
// --- Case 2: prestrain too long (positive) -> compression. Pinned both ends (stiff). ---
{
Model m; m.add_node(0,0); m.add_node(2,0);
double pre=1e-3; // too long
m.add_element(0,1,Section{E,A},MemberKind::Bar,pre);
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();
// fully restrained -> N = -EA*eps_free
near(r.axial_forces[0], -EA*pre, std::fabs(EA*pre)*1e-3, "C2 prestrain compression");
}
// --- Case 3: thermal heating dT>0 -> compression (restrained) ---
{
Model m; m.add_node(0,0); m.add_node(2,0);
double alpha=1.2e-5,dT=50;
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();
near(r.axial_forces[0], -EA*alpha*dT, std::fabs(EA*alpha*dT)*1e-3,"C3 thermal compression");
}
// --- Case 4: Cable in compression -> slack ---
{
Model m; 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_load(2,-1e4); // push node1 toward node0 -> compression
auto r=StaticSolver(m).solve();
chk(!r.active[0],"C4 cable slack");
near(r.axial_forces[0],0.0,1e-9,"C4 cable force 0");
near(r.utilization[0],0.0,1e-12,"C4 cable util 0");
}
// --- Case 4b: Cable in tension -> active ---
{
Model m; 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_load(2,1e4);
auto r=StaticSolver(m).solve();
chk(r.active[0],"C4b cable active tension");
near(r.axial_forces[0],1e4,1.0,"C4b cable tension=P");
}
// --- Case 5: Strut in tension -> separated ---
{
Model m; 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_load(2,1e4); // tension
auto r=StaticSolver(m).solve();
chk(!r.active[0],"C5 strut separated");
near(r.axial_forces[0],0.0,1e-9,"C5 strut force 0");
}
// --- Case 6: settlement loads structure with no applied force ---
{
Model m; m.add_node(0,0); m.add_node(2,0);
m.add_element(0,1,Section{E,A});
double K=1e9;
// node0 fixed stiff, node1: x-spring with settlement, y stiff
m.add_spring(0,1,0,1e16); m.add_spring(0,0,1,1e16); m.add_spring(1,0,1,1e16);
double s0=1e-3;
m.add_spring(1,1,0,K,s0); // settles +x by s0
auto r=StaticSolver(m).solve();
// 1D: node1 x. stiffness: bar k_b=EA/L, spring K. RHS = K*s0.
double kb=EA/2.0;
double u1 = K*s0/(kb+K);
near(r.displacements[2], u1, std::fabs(u1)*1e-6,"C6 settlement disp");
// axial = kb*u1 (tension, node1 pulled +x)
near(r.axial_forces[0], kb*u1, std::fabs(kb*u1)*1e-6,"C6 settlement axial");
// spring reaction = -K*(u1 - s0)
near(r.reactions[2], -K*(u1-s0), std::fabs(K*(u1-s0))*1e-6,"C6 spring reaction");
}
// --- Case 7: global equilibrium residual (Fext + reactions = internal at every node) ---
// Build a small 2-bar frame, solve, check sum of nodal equilibrium.
{
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,1e15); m.add_spring(0,0,1,1e15); m.add_spring(1,0,1,1e15);
m.add_load(4,-2e4); // node2 x
m.add_load(5,-3e4); // node2 y
auto r=StaticSolver(m).solve();
// nodal equilibrium: for each node, sum over members N*b_dir + reaction + load = 0
int nn=3;
std::vector<double> resid(2*nn,0.0);
for(int i=0;i<2*nn;++i) resid[i]= -r.reactions[i]; // start
for(auto&ld:m.loads()) resid[ld.dof]+= ld.value;
// member internal force on nodes = -N*b ... actually external eq: sum_members N*b = Fext+react
// We'll check: sum_members (N*b) - Fext - reactions == 0
std::vector<double> lhs(2*nn,0.0);
for(size_t e=0;e<m.num_elements();++e){
auto d=m.element_direction(e); auto&el=m.elements()[e];
double N=r.axial_forces[e];
double b[4]={-d.x,-d.y,d.x,d.y};
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) lhs[g[i]]+= N*b[i];
}
double maxr=0;
for(int i=0;i<2*nn;++i){ double rr=lhs[i]-r.reactions[i];
// also subtract applied loads
}
for(auto&ld:m.loads()) lhs[ld.dof]-=ld.value;
for(int i=0;i<2*nn;++i){ double rr=lhs[i]-r.reactions[i]; maxr=std::max(maxr,std::fabs(rr)); }
near(maxr,0.0,1e-3,"C7 global nodal equilibrium residual");
}
printf("\n%s (fails=%d)\n", fails?"SOME FAILED":"ALL PASSED", fails);
return fails?1:0;
}
EOF
g++ -std=c++17 -I/opt/truss2d/include /tmp/verify.cpp /opt/truss2d/src/*.cpp -o /tmp/verify 2>&1 | head -20 && /tmp/verifyExit code 127 /usr/bin/ld: /tmp/ccTXw4Vp.o: in function `main': main.cpp:(.text+0x0): multiple definition of `main'; /tmp/ccdZaP3O.o:verify.cpp:(.text+0x10b): first defined here collect2: error: ld returned 1 exit status /bin/bash: line 305: /tmp/verify: No such file or directory [error] tool reported failure
cd /opt/truss2d && g++ -std=c++17 -Iinclude /tmp/verify.cpp src/solver.cpp src/model.cpp src/dense_matrix.cpp -o /tmp/verify 2>&1 | head && /tmp/verify
Exit 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 1: 1746 Aborted /tmp/verify [error] tool reported failure
cd /opt/truss2d && cat > /tmp/v2.cpp <<'EOF'
#include <cmath>
#include <cstdio>
#include <vector>
#include "truss2d/model.hpp"
#include "truss2d/solver.hpp"
using namespace truss2d;
int fails=0;
void chk(bool c,const char*m){ if(!c){printf("FAIL: %s\n",m);++fails;} else printf("ok: %s\n",m);}
void near(double g,double w,double tol,const char*m){ if(std::fabs(g-w)>tol){printf("FAIL: %s got %.10g want %.10g\n",m,g,w);++fails;} else printf("ok: %s (%.6g)\n",m,g);}
int main(){
double E=2e11,A=1e-3,EA=E*A;
// Case 4: cable - weak x backup spring so slack case is determinate
{
Model m; m.add_node(0,0); m.add_node(2,0);
m.add_element(0,1,Section{E,A},MemberKind::Cable);
m.add_spring(0,1,0,1e15); m.add_spring(0,0,1,1e15); m.add_spring(1,0,1,1e15);
m.add_spring(1,1,0,1e3); // weak x backup
m.add_load(2,-1e4); // compression on cable
auto r=StaticSolver(m).solve();
chk(!r.active[0],"C4 cable slack under compression");
near(r.axial_forces[0],0,1e-9,"C4 force 0");
// node1 x carried entirely by weak spring: u=-1e4/1e3=-10
near(r.displacements[2],-10.0,1e-6,"C4 disp on weak spring");
}
// Case 4b tension
{
Model m; m.add_node(0,0); m.add_node(2,0);
m.add_element(0,1,Section{E,A},MemberKind::Cable);
m.add_spring(0,1,0,1e15); m.add_spring(0,0,1,1e15); m.add_spring(1,0,1,1e15);
m.add_spring(1,1,0,1e3);
m.add_load(2,1e4);
auto r=StaticSolver(m).solve();
chk(r.active[0],"C4b cable active");
// shares: kb=EA/2=1e8 huge vs 1e3 -> nearly all in cable
near(r.axial_forces[0],1e4,5.0,"C4b cable ~ P");
}
// Case 5 strut tension -> separated, weak backup
{
Model m; m.add_node(0,0); m.add_node(2,0);
m.add_element(0,1,Section{E,A},MemberKind::Strut);
m.add_spring(0,1,0,1e15); m.add_spring(0,0,1,1e15); m.add_spring(1,0,1,1e15);
m.add_spring(1,1,0,1e3);
m.add_load(2,1e4); // tension
auto r=StaticSolver(m).solve();
chk(!r.active[0],"C5 strut separated");
near(r.axial_forces[0],0,1e-9,"C5 force 0");
near(r.displacements[2],10.0,1e-6,"C5 disp on weak spring");
}
// Case 5b strut compression -> active
{
Model m; m.add_node(0,0); m.add_node(2,0);
m.add_element(0,1,Section{E,A},MemberKind::Strut);
m.add_spring(0,1,0,1e15); m.add_spring(0,0,1,1e15); m.add_spring(1,0,1,1e15);
m.add_spring(1,1,0,1e3);
m.add_load(2,-1e4);
auto r=StaticSolver(m).solve();
chk(r.active[0],"C5b strut active compression");
near(r.axial_forces[0],-1e4,5.0,"C5b strut ~ -P");
}
// Case 8: skew roller single spring at 45deg. Pin node0, node1 on 45-deg roller, load.
{
Model m; m.add_node(0,0); m.add_node(2,0);
m.add_element(0,1,Section{E,A});
m.add_spring(0,1,0,1e15); m.add_spring(0,0,1,1e15);
double inv=1.0/std::sqrt(2.0);
m.add_spring(1,inv,inv,1e15); // resists only along (1,1)/sqrt2
m.add_load(3,-1e4); // node1 y load
auto r=StaticSolver(m).solve();
// skew roller: no resistance perpendicular to (1,1). reaction must be along (1,1).
// reaction vector should be parallel to (1,1): Rx==Ry
near(r.reactions[2],r.reactions[3],std::fabs(r.reactions[2])*1e-6+1e-3,"C8 reaction parallel to spring dir (Rx=Ry)");
chk(std::fabs(r.reactions[2])>1.0,"C8 nonzero reaction");
}
// Case 9: geometric/initial-stress: pretensioned transverse cable net.
// Two cables along x pretensioned, node in middle, transverse y load.
// Without geometric stiffness, y-DOF has no elastic stiffness from axial cables (transverse).
{
Model m;
m.add_node(-1,0); m.add_node(1,0); m.add_node(0,0); // n2 center
// cables from center to anchors, pretensioned (manufactured short => tension): prestrain<0
double pre=-1e-3; // tension
m.add_element(2,0,Section{E,A},MemberKind::Cable,pre);
m.add_element(2,1,Section{E,A},MemberKind::Cable,pre);
// anchors fixed
m.add_spring(0,1,0,1e16); m.add_spring(0,0,1,1e16);
m.add_spring(1,1,0,1e16); m.add_spring(1,0,1,1e16);
// center: provide tiny x grounding to avoid x singularity; y resisted by geometric only
m.add_spring(2,1,0,1e2);
m.add_spring(2,0,1,1e-3); // negligible y spring just to be safe? keep tiny
double P=10.0; m.add_load(5,-P); // small transverse load at center y
auto r=StaticSolver(m).solve();
chk(r.active[0]&&r.active[1],"C9 both cables taut");
// pretension force magnitude ~ EA*|pre| = 2e8*1e-3=2e5
double T=EA*1e-3;
// transverse stiffness ~ sum 2 * T/L (L=1) = 2T ; plus tiny y spring 1e-3
double ky=2*T/1.0 + 1e-3;
double uy = -P/ky;
near(r.displacements[5],uy,std::fabs(uy)*1e-2,"C9 transverse disp from geometric stiffness");
chk(std::fabs(r.displacements[5])>0,"C9 nonzero transverse disp");
printf(" (T=%.6g, uy=%.6g, got=%.6g)\n",T,uy,r.displacements[5]);
}
printf("\n%s (fails=%d)\n", fails?"SOME FAILED":"ALL PASSED", fails);
return fails?1:0;
}
EOF
g++ -std=c++17 -Iinclude /tmp/v2.cpp src/solver.cpp src/model.cpp src/dense_matrix.cpp -o /tmp/v2 2>&1|head && /tmp/v2ok: C4 cable slack under compression ok: C4 force 0 (0) ok: C4 disp on weak spring (-10) ok: C4b cable active ok: C4b cable ~ P (9999.9) ok: C5 strut separated ok: C5 force 0 (0) ok: C5 disp on weak spring (10) ok: C5b strut active compression ok: C5b strut ~ -P (-9999.9) ok: C8 reaction parallel to spring dir (Rx=Ry) (9999.5) ok: C8 nonzero reaction ok: C9 both cables taut ok: C9 transverse disp from geometric stiffness (-2.5e-05) ok: C9 nonzero transverse disp (T=200000, uy=-2.5e-05, got=-2.5e-05) ALL PASSED (fails=0) [stdout] ok: C4 cable slack under compression ok: C4 force 0 (0) ok: C4 disp on weak spring (-10) ok: C4b cable active ok: C4b cable ~ P (9999.9) ok: C5 strut separated ok: C5 force 0 (0) ok: C5 disp on weak spring (10) ok: C5b strut active compression ok: C5b strut ~ -P (-9999.9) ok: C8 reaction parallel to spring dir (Rx=Ry) (9999.5) ok: C8 nonzero reaction ok: C9 both cables taut ok: C9 transverse disp from geometric stiffness (-2.5e-05) ok: C9 nonzero transverse disp (T=200000, uy=-2.5e-05, got=-2.5e-05) ALL PASSED (fails=0)
cat > /tmp/v3.cpp <<'EOF'
#include <cmath>
#include <cstdio>
#include <vector>
#include "truss2d/model.hpp"
#include "truss2d/solver.hpp"
using namespace truss2d;
int fails=0;
void near(double g,double w,double tol,const char*m){ if(std::fabs(g-w)>tol){printf("FAIL: %s got %.10g want %.10g\n",m,g,w);++fails;} else printf("ok: %s (%.6g)\n",m,g);}
int main(){
double E=2e11,A=1e-3,EA=E*A;
// C1 basic
{ Model m;m.add_node(0,0);m.add_node(2,0);m.add_element(0,1,Section{E,A});
m.add_spring(0,1,0,1e15);m.add_spring(0,0,1,1e15);m.add_spring(1,0,1,1e15);
m.add_load(2,1e4); auto r=StaticSolver(m).solve();
near(r.axial_forces[0],1e4,1.0,"C1 axial=P"); near(r.reactions[0],-1e4,1.0,"C1 Rx=-P");}
// C2 prestrain too long -> compression
{ Model m;m.add_node(0,0);m.add_node(2,0);double pre=1e-3;
m.add_element(0,1,Section{E,A},MemberKind::Bar,pre);
m.add_spring(0,1,0,1e16);m.add_spring(0,0,1,1e16);m.add_spring(1,1,0,1e16);m.add_spring(1,0,1,1e16);
auto r=StaticSolver(m).solve(); near(r.axial_forces[0],-EA*pre,std::fabs(EA*pre)*2e-3,"C2 too-long compression");}
// C3 thermal heat -> compression
{ Model m;m.add_node(0,0);m.add_node(2,0);double al=1.2e-5,dT=50;
m.add_element(0,1,Section{E,A},MemberKind::Bar,0,al,dT);
m.add_spring(0,1,0,1e16);m.add_spring(0,0,1,1e16);m.add_spring(1,1,0,1e16);m.add_spring(1,0,1,1e16);
auto r=StaticSolver(m).solve(); near(r.axial_forces[0],-EA*al*dT,std::fabs(EA*al*dT)*2e-3,"C3 heat compression");}
// C6 settlement
{ Model m;m.add_node(0,0);m.add_node(2,0);m.add_element(0,1,Section{E,A});
m.add_spring(0,1,0,1e16);m.add_spring(0,0,1,1e16);m.add_spring(1,0,1,1e16);
double K=1e9,s0=1e-3; m.add_spring(1,1,0,K,s0);
auto r=StaticSolver(m).solve(); double kb=EA/2,u1=K*s0/(kb+K);
near(r.displacements[2],u1,std::fabs(u1)*1e-6,"C6 settle disp");
near(r.axial_forces[0],kb*u1,std::fabs(kb*u1)*1e-6,"C6 settle axial");
near(r.reactions[2],-K*(u1-s0),std::fabs(K*(u1-s0))*1e-6,"C6 settle reaction");}
// C7 global equilibrium of a 3-bar truss
{ 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,1e15);m.add_spring(0,0,1,1e15);m.add_spring(1,0,1,1e15);
m.add_load(4,-2e4);m.add_load(5,-3e4); auto r=StaticSolver(m).solve();
std::vector<double> lhs(6,0.0);
for(size_t e=0;e<3;++e){auto d=m.element_direction(e);auto&el=m.elements()[e];double N=r.axial_forces[e];
double b[4]={-d.x,-d.y,d.x,d.y};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)lhs[g[i]]+=N*b[i];}
for(auto&ld:m.loads())lhs[ld.dof]-=ld.value;
double mx=0;for(int i=0;i<6;++i)mx=std::max(mx,std::fabs(lhs[i]-r.reactions[i]));
near(mx,0,1e-2,"C7 nodal equilibrium residual");}
printf("\n%s (fails=%d)\n",fails?"SOME FAILED":"ALL PASSED",fails); return fails?1:0;
}
EOF
g++ -std=c++17 -Iinclude /tmp/v3.cpp src/solver.cpp src/model.cpp src/dense_matrix.cpp -o /tmp/v3 2>&1|head && /tmp/v3Exit code 1 ok: C1 axial=P (10000) ok: C1 Rx=-P (-10000) ok: C2 too-long compression (-200000) ok: C3 heat compression (-120000) ok: C6 settle disp (0.000909091) ok: C6 settle axial (90909.1) ok: C6 settle reaction (90909.1) FAIL: C7 nodal equilibrium residual got 6.419094013 want 0 SOME FAILED (fails=1) [error] tool reported failure
cat > /tmp/v4.cpp <<'EOF'
#include <cmath>
#include <cstdio>
#include <vector>
#include "truss2d/model.hpp"
#include "truss2d/solver.hpp"
using namespace truss2d;
int main(){
double E=2e11,A=1e-3;
for(double Ks : {1e15, 1e11, 1e9}){
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,Ks);m.add_spring(0,0,1,Ks);m.add_spring(1,0,1,Ks);
m.add_load(4,-2e4);m.add_load(5,-3e4); auto r=StaticSolver(m).solve();
std::vector<double> lhs(6,0.0);
for(size_t e=0;e<3;++e){auto d=m.element_direction(e);auto&el=m.elements()[e];double N=r.axial_forces[e];
double b[4]={-d.x,-d.y,d.x,d.y};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)lhs[g[i]]+=N*b[i];}
for(auto&ld:m.loads())lhs[ld.dof]-=ld.value;
double mx=0,node2=0;for(int i=0;i<6;++i){double rr=std::fabs(lhs[i]-r.reactions[i]);mx=std::max(mx,rr); if(i>=4)node2=std::max(node2,rr);}
printf("Ks=%.0e max_resid=%.4g node2(free)_resid=%.4g\n",Ks,mx,node2);
}
}
EOF
g++ -std=c++17 -Iinclude /tmp/v4.cpp src/solver.cpp src/model.cpp src/dense_matrix.cpp -o /tmp/v4 && /tmp/v4Ks=1e+15 max_resid=6.419 node2(free)_resid=6.419 Ks=1e+11 max_resid=6.421 node2(free)_resid=6.421 Ks=1e+09 max_resid=6.653 node2(free)_resid=6.653 [stdout] Ks=1e+15 max_resid=6.419 node2(free)_resid=6.419 Ks=1e+11 max_resid=6.421 node2(free)_resid=6.421 Ks=1e+09 max_resid=6.653 node2(free)_resid=6.653
cat > /tmp/v5.cpp <<'EOF'
#include <cmath>
#include <cstdio>
#include <vector>
#include "truss2d/model.hpp"
#include "truss2d/solver.hpp"
using namespace truss2d;
int main(){
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_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,1e11);m.add_spring(0,0,1,1e11);m.add_spring(1,0,1,1e11);
m.add_load(4,-2e4);m.add_load(5,-3e4);
auto r=StaticSolver(m).solve();
auto&u=r.displacements;
std::vector<double> lhs(6,0.0);
for(size_t e=0;e<3;++e){auto d=m.element_direction(e);auto&el=m.elements()[e];double N=r.axial_forces[e];
double L=m.element_length(e); double c=d.x,s=d.y;
size_t g[4]={2*el.n1,2*el.n1+1,2*el.n2,2*el.n2+1};
double b[4]={-c,-s,c,s};
// axial internal:
for(int i=0;i<4;++i)lhs[g[i]]+=N*b[i];
// transverse geometric internal: K_g u = (N/L) a (a.u), a=(s,-c,-s,c)
double a[4]={s,-c,-s,c};
double au=0; for(int i=0;i<4;++i) au+=a[i]*u[g[i]];
for(int i=0;i<4;++i) lhs[g[i]] += (N/L)*a[i]*au;
}
for(auto&ld:m.loads())lhs[ld.dof]-=ld.value;
double mx=0;for(int i=0;i<6;++i)mx=std::max(mx,std::fabs(lhs[i]-r.reactions[i]));
printf("full equilibrium (axial+geometric) max residual = %.6g (loads ~3e4)\n",mx);
return mx<1e-4?0:1;
}
EOF
g++ -std=c++17 -Iinclude /tmp/v5.cpp src/solver.cpp src/model.cpp src/dense_matrix.cpp -o /tmp/v5 && /tmp/v5 && echo "EQUILIBRIUM VERIFIED"full equilibrium (axial+geometric) max residual = 3.63798e-12 (loads ~3e4) EQUILIBRIUM VERIFIED [stdout] full equilibrium (axial+geometric) max residual = 3.63798e-12 (loads ~3e4) EQUILIBRIUM VERIFIED
rm -rf build && cmake -S . -B build -DCMAKE_BUILD_TYPE=Release >/tmp/c.log 2>&1 && cmake --build build --parallel >/tmp/b.log 2>&1 && echo "=== BUILD OK ===" && cd build && ctest --output-on-failure 2>&1 | tail -8
=== 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 secNo 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.WOdWfvkpvX/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/model.cpp.o
[ 21%] Building CXX object CMakeFiles/truss2d.dir/src/io.cpp.o
[ 26%] Linking CXX static library libtruss2d.a
[ 26%] Built target truss2d
[ 31%] Building CXX object CMakeFiles/truss2d_cli.dir/src/main.cpp.o
[ 36%] Building CXX object tests/CMakeFiles/test_core.dir/test_core.cpp.o
[ 42%] Building CXX object tests/CMakeFiles/test_unilateral.dir/test_unilateral.cpp.o
[ 47%] Building CXX object tests/CMakeFiles/test_geometric.dir/test_geometric.cpp.o
[ 52%] Building CXX object tests/CMakeFiles/test_soak.dir/test_soak.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
[ 84%] Linking CXX executable test_soak
[ 84%] Linking CXX executable test_settlement
[ 89%] Linking CXX executable test_unilateral
[ 89%] Built target test_soak
[ 89%] Built target test_settlement
[ 94%] Linking CXX executable test_core
[ 94%] Built target test_unilateral
[ 94%] Built target test_core
[100%] Linking CXX executable test_geometric
[100%] Built target test_geometric
== Running hidden tests ==
Test project /tmp/tmp.WOdWfvkpvX/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_f13e993ee80a400c. Re-running the agent live requires EVAL_PLATFORM_ENABLE_OAUTH_SMOKE=1 and is non-deterministic.
Trial trial_f13e993ee80a400c · verifier authoritative; classifier explanatory.