SyncValsverifier → artifact → classifier → verdict
SyncVals · Trajectory

truss2d-solver

claude-code claude-opus-4-8 ✓ resolved GOOD_SUCCESS ↑ View task
Solved from the instruction alone, tests/ and solution/ were withheld from the agent's workspace and restored only for grading.
Reward = tests/test.sh exit code (0 → resolved); the classification below is post-hoc and cannot change it.
Classification , post-hoc; cannot change the reward
GOOD_SUCCESSLegitimate solve, implements the asked-for behavior; tests verify real functionality.
SubtypeLegitimate Solution
EvidenceTest output shows all 6 hidden tests passed: test_core, test_unilateral, test_geometric, test_soak, test_degenerate, test_settlement. Result.json confirms reward=1.0. Clean compilation with no build errors or warnings. Agent successfully implemented element_stiffness_global(), assemble_system(), prestrain_force(), and solve() functions meeting the complex structural solver specification.
Root causeThe agent (Claude Opus 4.8) correctly implemented a complex finite element structural solver handling unilateral members (cables/struts), prestrain, thermal strain, spring supports, and initial-stress stiffness effects as specified in the detailed instruction.md.
RecommendationN/A - task is fine
Trajectory
Tool-by-tool agent trajectory
26 tool calls · 4 tool types · 44 steps
# Ticket: Implement the static solver for the `truss2d` network engine ## Context `truss2d` is a compact C++17 structural-analysis engine for planar networks of two-force members on elastic spring supports. The framework around the numerical core is already implemented: the data model and validation, dense linear algebra, text parser, report writer, demo CLI, and a smoke test. The static solver itself is intentionally stubbed. It returns zero-valued results, so the engine does not carry load or satisfy equilibrium. The project is at **`/opt/truss2d`** in the build image. ## Your Task Implement the four member functions in **`/opt/truss2d/src/solver.cpp`** (declared in `include/truss2d/solver.hpp`): - `DenseMatrix StaticSolver::element_stiffness_global(std::size_t e) const` - `DenseMatrix StaticSolver::assemble_system(const std::vector<bool>& active) const` - `std::vector<double> StaticSolver::prestrain_force(const std::vector<bool>& active) const` - `SolveResult StaticSolver::solve() const` Do not change public headers or signatures. You should only need to edit `src/solver.cpp`, using the existing `Model`, `DenseMatrix`, and `solve_spd` support code. This is not just a classical pin-jointed linear truss. The solver must handle unilateral members, installed strain, thermal strain, finite spring supports, and the small-displacement initial-stress effect of taut members. ## Governing Model Each node has two translational degrees of freedom, `x` then `y`. Each element connects two nodes and acts as a two-force member in its own axis. The solver uses small-displacement kinematics: axial extension is the first-order change of end separation projected onto the member axis. ### Member Response Each member has a linear-elastic axial response based on its material, area, and undeformed length. The `kind` controls which signs of total axial force it may carry: - `Bar`: carries either tension or compression and always participates. - `Cable`: tension-only. If the consistent state would place it in compression, it goes slack and contributes no stiffness or installed load. - `Strut`: compression-only. If the consistent state would place it in tension, it separates and contributes no stiffness or installed load. The sign convention for `SolveResult::axial_forces` is positive in tension and negative in compression. Dropped unilateral members must report exactly zero force and zero utilization. ### Installed Strain `Element::prestrain` is a mechanical lack-of-fit strain. `Element::alpha` and `Element::dT` describe a thermal strain contribution. Combine the mechanical and thermal installed strains consistently to first order before computing the installed axial force. A member manufactured too long for its joints pushes its ends apart and is initially compressive. A member manufactured too short pulls its ends together and is initially tensile. The same physical convention applies to thermal expansion and contraction. ### Initial-Stress Stiffness An active member carrying axial force also changes the tangent stiffness for relative motion transverse to its axis. Tension stiffens the transverse mode; compression softens it. This initial-stress contribution: - uses the member's current total axial force, not just the installed component; - acts only for participating members; - does not alter the recovered axial force itself; - vanishes for rigid translation of both ends and for relative motion along the member axis. The transverse stiffness of a taut member depends on the force it carries, and the force depends on the displacement. A single linear solve made from an arbitrary assumed member force is not, by itself, an equilibrium. Find the state in which the assumed member forces, the displacements they produce, and the forces recovered from those displacements all coincide. The graded models are well posed and have a unique consistent state. The matrix helpers report the ordinary elastic behavior only: the single-member helper reports one member's ordinary elastic stiffness, and the assembly helper adds the ordinary elastic stiffness of participating members together with the spring-support stiffness. The force-dependent transverse effect belongs only to the equilibrium solved by `solve()`, not to those helper outputs. ### Spring Supports There are no exact fixed supports. Each support is a grounded linear spring at a node along a unit direction stored by the model. It resists only the displacement component along that direction and produces no reaction perpendicular to it. Two nonparallel springs on one node act like a two-dimensional elastic support; one spring is a skew roller with finite compliance. The support reaction reported by `solve()` is the force applied by the spring to the structure, so it opposes the spring extension. A spring support may also **settle**: its grounded end is held at a prescribed offset (the `settlement` value, measured along the spring's direction). The spring then stores its stretch RELATIVE to that offset , i.e. the elongation is the node's displacement along the direction minus the prescribed offset , so a nonzero settlement loads the structure even with no applied forces, and the reaction depends on that relative stretch. A settlement of zero is an ordinary grounded spring. ### Equilibrium And Active Set At convergence, every node is in static equilibrium under applied loads, support reactions, the axial forces of active members, and the transverse initial-stress forces implied by those active members. The active set and the displacement field must agree with the unilateral rules above. A cable or strut is judged by its total axial force after installed strain and elastic extension are both included. If a governing system is singular or not positive definite, let the exception from the provided linear solver propagate. `solve()` must call `Model::validate()` before assembling the analysis. ## Outputs `SolveResult` must populate: - `displacements`: one value per global DOF. - `reactions`: one value per global DOF, zero for unsupported DOFs. - `axial_forces`: one value per element, positive in tension. - `utilization`: one signed value per element. Use an allowable stress of `2.5e8 Pa`; capacity scales with area. Inactive unilateral members report zero. - `active`: one flag per element. All vectors must have the lengths implied by the model. ## Build And Test ```bash cmake -S /opt/truss2d -B /opt/truss2d/build -DCMAKE_BUILD_TYPE=Release cmake --build /opt/truss2d/build --parallel cd /opt/truss2d/build && ctest --output-on-failure ``` The shipped test is only a smoke test. The grader restores a clean project snapshot, copies in your `src/solver.cpp`, and runs a hidden correctness suite against independent analytic and numerical references.
I'll start by exploring the project structure and understanding the existing code.

/opt/truss2d/include/truss2d/solver.hpp

contents
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

contents
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
Let me read the key support files to understand the data model and helpers.

/opt/truss2d/include/truss2d/model.hpp

contents
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

contents
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

contents
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
Let me look at the model implementation, instruction.md, and existing tests.

/opt/truss2d/src/model.cpp

contents
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

Read instruction.md and list project files

shell
cat /opt/truss2d/instruction.md 2>/dev/null; echo "---FILES---"; ls -R /opt/truss2d --ignore=build
output
---FILES---
/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]
---FILES---
/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
Let me read the tests, the report writer, and the demo to understand usage conventions.

/opt/truss2d/tests/test_smoke.cpp

contents
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

contents
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/dense_matrix.cpp

contents
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
Let me check io.cpp and main.cpp for output conventions like utilization sign.

/opt/truss2d/src/io.cpp

contents
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

contents
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
I now understand the full model. Let me verify there's no `instruction.md` with extra detail and check the example file, then implement.

Check for instruction.md and example

shell
ls -la /opt/truss2d/*.md 2>/dev/null; echo "---"; cat /opt/truss2d/examples/two_bar.truss; echo "---CMake tests---"; cat /opt/truss2d/tests/CMakeLists.txt
output
---
# Example network for the truss2d demo CLI.
#
# A free node (2) is suspended below two top anchors (0,1) by two CABLES and
# propped against a ground node (3) by a BAR. A sideways pull at node 2 slackens
# one cable; the other cable plus the bar carry the load.
#
#   node <x> <y>
node -3.0  4.0     # 0 top-left anchor
node  3.0  4.0     # 1 top-right anchor
node  0.0  0.0     # 2 free node (loaded)
node  0.0 -4.0     # 3 ground

# elem <n1> <n2> <E> <A> [kind: 0=bar 1=cable 2=strut] [prestrain]
elem 2 0 2.0e11 1.0e-3 1        # cable up-left
elem 2 1 2.0e11 1.0e-3 1        # cable up-right
elem 2 3 2.0e11 1.0e-3 0        # bar down to ground

# spring <node> <dx> <dy> <kappa>   (grounded spring along (dx,dy), N/m)
spring 0 1 0 5.0e10
spring 0 0 1 5.0e10
spring 1 1 0 5.0e10
spring 1 0 1 5.0e10
spring 3 1 0 5.0e10
spring 3 0 1 5.0e10

# load <node> <comp> <value>   (comp: 0 = x, 1 = y)
load 2 0 4.0e4     # push node 2 to the right
---CMake tests---
# Smoke test for the shipped engine. Each test file is a standalone executable
# using the tiny header-only harness in test_util.hpp. CTest runs them and
# treats a non-zero exit as failure.
#
# NOTE: this directory ships only a trivial smoke test. The full correctness /
# invariant / discriminator suite is supplied by the grader at evaluation time
# (it adds its own test sources and registers them with CTest).

add_executable(test_smoke test_smoke.cpp)
target_link_libraries(test_smoke PRIVATE truss2d)
target_include_directories(test_smoke PRIVATE ${CMAKE_CURRENT_SOURCE_DIR})
add_test(NAME smoke COMMAND test_smoke)

[stdout]
---
# Example network for the truss2d demo CLI.
#
# A free node (2) is suspended below two top anchors (0,1) by two CABLES and
# propped against a ground node (3) by a BAR. A sideways pull at node 2 slackens
# one cable; the other cable plus the bar carry the load.
#
#   node <x> <y>
node -3.0  4.0     # 0 top-left anchor
node  3.0  4.0     # 1 top-right anchor
node  0.0  0.0     # 2 free node (loaded)
node  0.0 -4.0     # 3 ground

# elem <n1> <n2> <E> <A> [kind: 0=bar 1=cable 2=strut] [prestrain]
elem 2 0 2.0e11 1.0e-3 1        # cable up-left
elem 2 1 2.0e11 1.0e-3 1        # cable up-right
elem 2 3 2.0e11 1.0e-3 0        # bar down to ground

# spring <node> <dx> <dy> <kappa>   (grounded spring along (dx,dy), N/m)
spring 0 1 0 5.0e10
spring 0 0 1 5.0e10
spring 1 1 0 5.0e10
spring 1 0 1 5.0e10
spring 3 1 0 5.0e10
spring 3 0 1 5.0e10

# load <node> <comp> <value>   (comp: 0 = x, 1 = y)
load 2 0 4.0e4     # push node 2 to the right
---CMake tests---
# Smoke test for the shipped engine. Each test file is a standalone executable
# using the tiny header-only harness in test_util.hpp. CTest runs them and
# treats a non-zero exit as failure.
#
# NOTE: this directory ships only a trivial smoke test. The full correctness /
# invariant / discriminator suite is supplied by the grader at evaluation time
# (it adds its own test sources and registers them with CTest).

add_executable(test_smoke test_smoke.cpp)
target_link_libraries(test_smoke PRIVATE truss2d)
target_include_directories(test_smoke PRIVATE ${CMAKE_CURRENT_SOURCE_DIR})
add_test(NAME smoke COMMAND test_smoke)
Now I have a complete understanding. Let me write the implementation.
write · /opt/truss2d/src/solver.cpp
#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 on spring supports.
//
// Conventions used throughout:
//   * Each node owns DOFs {2*node (x), 2*node + 1 (y)}.
//   * For an element from n1 to n2 with direction cosines (c, s), the axial
//     extension is the first-order change of end separation projected on the
//     member axis: delta = B . u, with B = [-c, -s, c, s] over the element's
//     four DOFs {n1.x, n1.y, n2.x, n2.y}.
//   * The recovered (total) axial force is N = (EA/L) * delta - EA * eps0,
//     with eps0 the installed (mechanical + thermal) stress-free strain. N is
//     positive in tension. A member manufactured too long (eps0 > 0) is
//     compressive at zero displacement; too short (eps0 < 0) is tensile.
//   * Elastic stiffness K0 = (EA/L) B^T B (axial only).
//   * Geometric (initial-stress) stiffness K_g = (N/L) Bn^T Bn, where
//     Bn = [s, -c, -s, c] projects relative motion onto the member normal
//     (-s, c). Tension (N > 0) stiffens the transverse mode; compression
//     softens it. K_g vanishes for rigid translation and for axial relative
//     motion, and never alters the recovered axial force itself.
// ============================================================================

namespace {

constexpr double kAllowableStress = 2.5e8; // [Pa]

// Per-element kinematic/material cache.
struct ElementInfo {
    double c{0.0};    // cos of axis
    double s{0.0};    // sin of axis
    double L{0.0};    // undeformed length
    double EA{0.0};   // axial rigidity
    double k{0.0};    // EA / L
    double eps0{0.0}; // installed stress-free strain (mechanical + thermal)
    std::array<std::size_t, 4> dof{}; // {n1.x, n1.y, n2.x, n2.y}
};

} // namespace

DenseMatrix StaticSolver::element_stiffness_global(std::size_t e) const {
    const Vec2 d = model_.element_direction(e);
    const double L = model_.element_length(e);
    const Element& el = model_.elements()[e];
    const double k = el.section.E * el.section.A / L;

    const std::array<double, 4> B{-d.x, -d.y, d.x, d.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 * 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);

    // Ordinary elastic stiffness of participating members.
    const std::size_t ne = model_.num_elements();
    for (std::size_t e = 0; e < ne; ++e) {
        if (e < active.size() && !active[e]) {
            continue;
        }
        const Element& el = model_.elements()[e];
        const std::array<std::size_t, 4> dof{
            2 * el.n1, 2 * el.n1 + 1, 2 * el.n2, 2 * el.n2 + 1};
        const DenseMatrix Ke = element_stiffness_global(e);
        for (std::size_t i = 0; i < 4; ++i) {
            for (std::size_t j = 0; j < 4; ++j) {
                K(dof[i], dof[j]) += Ke(i, j);
            }
        }
    }

    // Spring-support stiffness: kappa * d d^T at the supported node.
    for (const SpringSupport& sp : model_.springs()) {
        const std::size_t bx = 2 * sp.node;
        const std::size_t by = 2 * sp.node + 1;
        const double kxx = sp.kappa * sp.dir.x * sp.dir.x;
        const double kxy = sp.kappa * sp.dir.x * sp.dir.y;
        const double kyy = sp.kappa * sp.dir.y * sp.dir.y;
        K(bx, bx) += kxx;
        K(bx, by) += kxy;
        K(by, bx) += kxy;
        K(by, by) += kyy;
    }

    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 std::size_t ne = model_.num_elements();

    for (std::size_t e = 0; e < ne; ++e) {
        if (e < active.size() && !active[e]) {
            continue;
        }
        const Element& el = model_.elements()[e];
        const Vec2 d = model_.element_direction(e);
        const double eps0 = el.prestrain + el.alpha * el.dT;
        // Installed nodal force g = EA * eps0 applied along B^T: a member that
        // is too long (eps0 > 0) pushes its ends apart.
        const double g = el.section.E * el.section.A * eps0;
        f[2 * el.n1] += -g * d.x;
        f[2 * el.n1 + 1] += -g * d.y;
        f[2 * el.n2] += g * d.x;
        f[2 * el.n2 + 1] += g * d.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();

    // Cache per-element kinematics/material data.
    std::vector<ElementInfo> info(ne);
    for (std::size_t e = 0; e < ne; ++e) {
        const Element& el = model_.elements()[e];
        const Vec2 d = model_.element_direction(e);
        const double L = model_.element_length(e);
        ElementInfo& fo = info[e];
        fo.c = d.x;
        fo.s = d.y;
        fo.L = L;
        fo.EA = el.section.E * el.section.A;
        fo.k = fo.EA / L;
        fo.eps0 = el.prestrain + el.alpha * el.dT;
        fo.dof = {2 * el.n1, 2 * el.n1 + 1, 2 * el.n2, 2 * el.n2 + 1};
    }

    // Constant right-hand side parts that do not depend on the active set:
    // applied loads and spring settlement.
    std::vector<double> f_fixed(n, 0.0);
    for (const Load& ld : model_.loads()) {
        f_fixed[ld.dof] += ld.value;
    }
    for (const SpringSupport& sp : model_.springs()) {
        // Settlement moves the spring's grounded end; the stored stretch is
        // (u.dir - settlement), so a nonzero settlement loads the structure:
        // RHS contribution +kappa * settlement * dir.
        const double g = sp.kappa * sp.settlement;
        f_fixed[2 * sp.node] += g * sp.dir.x;
        f_fixed[2 * sp.node + 1] += g * sp.dir.y;
    }

    // Recover the axial force of element e from a displacement field.
    auto axial_force = [&](std::size_t e, const std::vector<double>& u) {
        const ElementInfo& fo = info[e];
        const double delta = -fo.c * u[fo.dof[0]] - fo.s * u[fo.dof[1]] +
                             fo.c * u[fo.dof[2]] + fo.s * u[fo.dof[3]];
        return fo.k * delta - fo.EA * fo.eps0;
    };

    std::vector<bool> active(ne, true);
    std::vector<double> N(ne, 0.0); // assumed member forces (geometric stiffness)
    std::vector<double> u(n, 0.0);

    const std::size_t max_outer = 4 * ne + 20;
    for (std::size_t outer = 0; outer < max_outer; ++outer) {
        // --- Inner fixed point: find forces consistent with the displacements
        //     they produce, for the current active set. ---
        const std::size_t max_inner = 200;
        for (std::size_t e = 0; e < ne; ++e) {
            if (!active[e]) N[e] = 0.0;
        }
        for (std::size_t inner = 0; inner < max_inner; ++inner) {
            DenseMatrix K = assemble_system(active);
            // Add the geometric (initial-stress) stiffness of active members
            // using the current assumed axial force.
            for (std::size_t e = 0; e < ne; ++e) {
                if (!active[e]) continue;
                const ElementInfo& fo = info[e];
                const double f = N[e] / fo.L;
                const std::array<double, 4> Bn{fo.s, -fo.c, -fo.s, fo.c};
                for (std::size_t i = 0; i < 4; ++i) {
                    for (std::size_t j = 0; j < 4; ++j) {
                        K(fo.dof[i], fo.dof[j]) += f * Bn[i] * Bn[j];
                    }
                }
            }

            std::vector<double> rhs = f_fixed;
            const std::vector<double> fp = prestrain_force(active);
            for (std::size_t i = 0; i < n; ++i) rhs[i] += fp[i];

            u = solve_spd(K, rhs);

            // Recover forces; measure the change from the assumed forces.
            double max_change = 0.0;
            double max_mag = 0.0;
            std::vector<double> N_new(ne, 0.0);
            for (std::size_t e = 0; e < ne; ++e) {
                if (!active[e]) continue;
                const double nf = axial_force(e, u);
                N_new[e] = nf;
                max_change = std::max(max_change, std::fabs(nf - N[e]));
                max_mag = std::max(max_mag, std::fabs(nf));
            }
            N = N_new;

            if (max_change <= 1e-11 * (max_mag + 1.0)) {
                break;
            }
        }

        // --- Update the active set from the trial forces of EVERY member,
        //     evaluated on the current displacement field. ---
        std::vector<bool> next = active;
        bool changed = false;
        for (std::size_t e = 0; e < ne; ++e) {
            const MemberKind kind = model_.elements()[e].kind;
            bool want = true;
            if (kind == MemberKind::Bar) {
                want = true;
            } else {
                const double trial = axial_force(e, u);
                if (kind == MemberKind::Cable) {
                    want = trial > 0.0; // tension-only
                } else {                // Strut: compression-only
                    want = trial < 0.0;
                }
            }
            if (want != active[e]) {
                next[e] = want;
                changed = true;
            }
        }

        if (!changed) {
            break;
        }
        active = next;
    }

    // --- Assemble the result from the converged state. ---
    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; // slack cable / separated strut: exactly zero
        }
        const double nf = axial_force(e, u);
        result.axial_forces[e] = nf;
        const double capacity = kAllowableStress * model_.elements()[e].section.A;
        result.utilization[e] = nf / capacity;
    }

    // Support reactions: the force applied by each spring to the structure,
    // opposing its (settlement-relative) extension.
    for (const SpringSupport& 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 force = -sp.kappa * ext;
        result.reactions[2 * sp.node] += force * sp.dir.x;
        result.reactions[2 * sp.node + 1] += force * sp.dir.y;
    }

    return result;
}

} // namespace truss2d
I need to add `<algorithm>` for `std::max`. Let me fix that and build.

/opt/truss2d/src/solver.cpp

− removed
#include <array>
#include <cmath>
+ added
#include <algorithm>
#include <array>
#include <cmath>

Build and run tests

shell
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
output
[ 11%] Building CXX object CMakeFiles/truss2d.dir/src/dense_matrix.cpp.o
[ 22%] Building CXX object CMakeFiles/truss2d.dir/src/io.cpp.o
[ 33%] Building CXX object CMakeFiles/truss2d.dir/src/model.cpp.o
[ 44%] Building CXX object CMakeFiles/truss2d.dir/src/solver.cpp.o
[ 55%] Linking CXX static library libtruss2d.a
[ 55%] Built target truss2d
[ 77%] 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/io.cpp.o
[ 33%] Building CXX object CMakeFiles/truss2d.dir/src/model.cpp.o
[ 44%] Building CXX object CMakeFiles/truss2d.dir/src/solver.cpp.o
[ 55%] Linking CXX static library libtruss2d.a
[ 55%] Built target truss2d
[ 77%] 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
Builds and the smoke test passes. Let me verify physical correctness against the example and check a few analytic cases with a scratch test.

Run demo CLI on example

shell
./truss2d_cli ../examples/two_bar.truss
output
# Truss2D static analysis report
nodes=4 elements=3 dofs=8

[displacements] (node, ux, uy)
0 7.988854e-07 -1.063410e-06
1 0.000000e+00 0.000000e+00
2 4.192473e-03 1.064474e-03
3 1.114579e-09 1.063410e-06

[reactions] (node, Rx, Ry)
0 -3.994427e+04 5.317052e+04
1 0.000000e+00 0.000000e+00
2 0.000000e+00 0.000000e+00
3 -5.572897e+01 -5.317052e+04

[axial_forces] (element, kind, active, N [+tension], utilization)
0 cable 1 6.650298e+04 2.660119e-01
1 cable 0 0.000000e+00 0.000000e+00
2 bar 1 5.317052e+04 2.126821e-01

[stdout]
# Truss2D static analysis report
nodes=4 elements=3 dofs=8

[displacements] (node, ux, uy)
0 7.988854e-07 -1.063410e-06
1 0.000000e+00 0.000000e+00
2 4.192473e-03 1.064474e-03
3 1.114579e-09 1.063410e-06

[reactions] (node, Rx, Ry)
0 -3.994427e+04 5.317052e+04
1 0.000000e+00 0.000000e+00
2 0.000000e+00 0.000000e+00
3 -5.572897e+01 -5.317052e+04

[axial_forces] (element, kind, active, N [+tension], utilization)
0 cable 1 6.650298e+04 2.660119e-01
1 cable 0 0.000000e+00 0.000000e+00
2 bar 1 5.317052e+04 2.126821e-01
The example behaves correctly (one cable slackens, reactions balance). Let me write a standalone verification harness for analytic cases: prestrain/thermal sign, settlement, geometric transverse stiffness, and unilateral logic.

Compile and run analytic verification

shell
cat > /tmp/verify.cpp <<'EOF'
#include <cstdio>
#include <cmath>
#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*(1+std::fabs(want));
  printf("[%s] %-32s got=%.6e want=%.6e\n", ok?"PASS":"FAIL",name,got,want);
  if(!ok) fails++;
}

int main(){
  // ---- Case 1: single horizontal bar, axial load. EA/L stiffness. ----
  // node0 fully sprung (stiff), node1 sprung in y only, load +P in x.
  {
    Model m; double E=2e11,A=1e-3,L=2.0,P=1e4;
    m.add_node(0,0); m.add_node(L,0);
    m.add_element(0,1,Section{E,A});
    double big=1e15;
    m.add_spring(0,1,0,big); m.add_spring(0,0,1,big); m.add_spring(1,0,1,big);
    m.add_load(2,P);
    auto r=StaticSolver(m).solve();
    double u=r.displacements[2];
    chk("bar axial disp", u, P*L/(E*A), 1e-6);
    chk("bar axial force", r.axial_forces[0], P, 1e-6); // tension
    chk("bar reaction x@0", r.reactions[0], -P, 1e-6);
  }

  // ---- Case 2: thermal heating of a bar restrained both ends -> compression ----
  {
    Model m; double E=2e11,A=1e-3,L=2.0,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 big=1e15;
    m.add_spring(0,1,0,big); m.add_spring(0,0,1,big);
    m.add_spring(1,1,0,big); m.add_spring(1,0,1,big);
    auto r=StaticSolver(m).solve();
    chk("thermal force (compress)", r.axial_forces[0], -E*A*alpha*dT, 1e-6);
  }

  // ---- Case 3: prestrain too-long (eps0>0) restrained -> compression ----
  {
    Model m; double E=2e11,A=1e-3,L=2.0,pre=1e-4;
    m.add_node(0,0); m.add_node(L,0);
    m.add_element(0,1,Section{E,A},MemberKind::Bar,pre);
    double big=1e15;
    m.add_spring(0,1,0,big); m.add_spring(0,0,1,big);
    m.add_spring(1,1,0,big); m.add_spring(1,0,1,big);
    auto r=StaticSolver(m).solve();
    chk("prestrain too-long compress", r.axial_forces[0], -E*A*pre, 1e-6);
  }

  // ---- Case 4: settlement. single x-spring at node1, member to node0(fixed).
  // settlement s along +x: spring grounded end moves +s. With member EA/L and
  // spring kappa in series. node1 free in x. ----
  {
    Model m; double E=2e11,A=1e-3,L=2.0,kappa=1e7,s=1e-3;
    m.add_node(0,0); m.add_node(L,0);
    m.add_element(0,1,Section{E,A});
    double big=1e15;
    m.add_spring(0,1,0,big); m.add_spring(0,0,1,big); m.add_spring(1,0,1,big);
    m.add_spring(1,1,0,kappa,s); // settling x-spring at node1
    auto r=StaticSolver(m).solve();
    double kbar=E*A/L;
    // series: u1 = kappa*s/(kbar+kappa)
    double u1=kappa*s/(kbar+kappa);
    chk("settlement disp", r.displacements[2], u1, 1e-6);
    double N=kbar*u1; // bar tension pulled by spring
    chk("settlement bar force", r.axial_forces[0], N, 1e-6);
    // spring reaction = -kappa*(u1 - s)
    chk("settlement reaction", r.reactions[2], -kappa*(u1-s), 1e-6);
  }

  // ---- Case 5: cable goes slack under compression-inducing load ----
  {
    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 big=1e15;
    m.add_spring(0,1,0,big); m.add_spring(0,0,1,big); m.add_spring(1,0,1,big);
    m.add_load(2,-1e4); // push node1 toward node0 -> compression -> slack
    auto r=StaticSolver(m).solve();
    chk("cable slack force", r.axial_forces[0], 0.0, 1e-9);
    printf("    cable active=%d (want 0)\n", (int)r.active[0]);
    if(r.active[0]) fails++;
  }

  // ---- Case 6: geometric transverse stiffness of a pretensioned cable ----
  // node0(0,0) fixed, node1(L,0) x-fixed stiff, free in y, transverse load P.
  // pretension via prestrain too-short (eps0<0 -> tension T0). Transverse
  // stiffness = T0/L, so v = P/(T0/L) = P*L/T0.
  {
    Model m; double E=2e11,A=1e-3,L=2.0,pre=-5e-4,P=50.0;
    m.add_node(0,0); m.add_node(L,0);
    m.add_element(0,1,Section{E,A},MemberKind::Bar,pre);
    double big=1e16;
    m.add_spring(0,1,0,big); m.add_spring(0,0,1,big); m.add_spring(1,1,0,big);
    m.add_load(3,P); // y at node1
    auto r=StaticSolver(m).solve();
    double T0=-E*A*pre; // tension
    chk("pretension force", r.axial_forces[0], T0, 1e-4);
    chk("transverse disp (geom)", r.displacements[3], P*L/T0, 2e-3);
  }

  printf("\n%s (%d failures)\n", fails? "FAILURES":"ALL PASS", 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 && /tmp/verify
output
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 223:  1760 Aborted                 /tmp/verify

[error] tool reported failure
Case 6 exposes a real bug: seeding the assumed force at **zero** makes the first solve singular for a transversely-loaded pretensioned member (its transverse stiffness comes entirely from the geometric term). I should seed the assumed force with the **installed force** `-EA·eps0` so pretensioned members have transverse stiffness from the first solve.

/opt/truss2d/src/solver.cpp

− removed
    std::vector<bool> active(ne, true);
    std::vector<double> N(ne, 0.0); // assumed member forces (geometric stiffness)
    std::vector<double> u(n, 0.0);
+ added
    std::vector<bool> active(ne, true);
    // Assumed member forces driving the geometric stiffness. Seed each member
    // with its installed force (-EA*eps0): a pretensioned member then carries
    // transverse stiffness from the very first solve, which a zero seed would
    // miss (leaving a transversely-loaded taut member singular).
    std::vector<double> N(ne, 0.0);
    for (std::size_t e = 0; e < ne; ++e) {
        N[e] = -info[e].EA * info[e].eps0;
    }
    std::vector<double> u(n, 0.0);
Now reseed members that transition from inactive to active (their stored force was zeroed), so they also get a proper geometric seed:

/opt/truss2d/src/solver.cpp

− removed
            if (want != active[e]) {
                next[e] = want;
                changed = true;
            }
+ added
            if (want != active[e]) {
                next[e] = want;
                changed = true;
                if (want) {
                    // Reactivated member: reseed its assumed force with the
                    // installed force so it has transverse stiffness on the
                    // next solve.
                    N[e] = -info[e].EA * info[e].eps0;
                }
            }

Rebuild and rerun verification

shell
cmake --build /opt/truss2d/build --parallel 2>&1 | tail -3 && g++ -std=c++17 -I/opt/truss2d/include /tmp/verify.cpp /opt/truss2d/build/libtruss2d.a -o /tmp/verify && /tmp/verify
output
Exit code 134
[ 77%] Linking CXX executable truss2d_cli
[ 88%] Built target truss2d_cli
[100%] Built target test_smoke
terminate called after throwing an instance of 'std::runtime_error'
  what():  solve_spd: matrix is singular or not positive-definite
/bin/bash: line 1:  1817 Aborted                 /tmp/verify

[error] tool reported failure
Case 5 is actually an ill-posed *test* (a slack cable leaves node1 unsupported in x → genuinely singular, which is correct behavior). Let me wrap each case in try/catch to see all results and fix that test.

Locate case 5

shell
cd /tmp && python3 - <<'PY'
import re
s=open('verify.cpp').read()
# wrap each block { ... } in main with try/catch by replacing the case openers
PY
# simpler: edit case5 to add x-spring at node1 and catch-all in main
sed -i 's|    m.add_spring(0,1,0,big); m.add_spring(0,0,1,big); m.add_spring(1,0,1,big);\n    m.add_load(2,-1e4);|XX|' verify.cpp
grep -n "push node1 toward node0" verify.cpp
output
83:    m.add_load(2,-1e4); // push node1 toward node0 -> compression -> slack
Shell cwd was reset to /opt/truss2d

[stdout]
83:    m.add_load(2,-1e4); // push node1 toward node0 -> compression -> slack
[stderr]

Shell cwd was reset to /opt/truss2d

Fix case 5 and rerun verification

shell
cd /tmp && python3 - <<'PY'
s=open('verify.cpp').read()
s=s.replace(
"""    m.add_spring(0,1,0,big); m.add_spring(0,0,1,big); m.add_spring(1,0,1,big);
    m.add_load(2,-1e4); // push node1 toward node0 -> compression -> slack""",
"""    double ksm=1e6;
    m.add_spring(0,1,0,big); m.add_spring(0,0,1,big); m.add_spring(1,0,1,big);
    m.add_spring(1,1,0,ksm); // finite x support at node1
    m.add_load(2,-1e4); // push node1 toward node0 -> compression -> slack""")
# add per-case try/catch is messy; instead wrap whole main body. Just trust now.
open('verify.cpp','w').write(s)
print("ok")
PY
g++ -std=c++17 -I/opt/truss2d/include /tmp/verify.cpp /opt/truss2d/build/libtruss2d.a -o /tmp/verify && /tmp/verify
output
ok
[PASS] bar axial disp                   got=1.000000e-04 want=1.000000e-04
[PASS] bar axial force                  got=1.000000e+04 want=1.000000e+04
[PASS] bar reaction x@0                 got=-1.000000e+04 want=-1.000000e+04
[PASS] thermal force (compress)         got=-1.200000e+05 want=-1.200000e+05
[PASS] prestrain too-long compress      got=-2.000000e+04 want=-2.000000e+04
[PASS] settlement disp                  got=9.090910e-05 want=9.090909e-05
[PASS] settlement bar force             got=9.090909e+03 want=9.090909e+03
[PASS] settlement reaction              got=9.090909e+03 want=9.090909e+03
[PASS] cable slack force                got=0.000000e+00 want=0.000000e+00
    cable active=0 (want 0)
[PASS] pretension force                 got=1.000000e+05 want=1.000000e+05
[PASS] transverse disp (geom)           got=1.000000e-03 want=1.000000e-03

ALL PASS (0 failures)
Shell cwd was reset to /opt/truss2d

[stdout]
ok
[PASS] bar axial disp                   got=1.000000e-04 want=1.000000e-04
[PASS] bar axial force                  got=1.000000e+04 want=1.000000e+04
[PASS] bar reaction x@0                 got=-1.000000e+04 want=-1.000000e+04
[PASS] thermal force (compress)         got=-1.200000e+05 want=-1.200000e+05
[PASS] prestrain too-long compress      got=-2.000000e+04 want=-2.000000e+04
[PASS] settlement disp                  got=9.090910e-05 want=9.090909e-05
[PASS] settlement bar force             got=9.090909e+03 want=9.090909e+03
[PASS] settlement reaction              got=9.090909e+03 want=9.090909e+03
[PASS] cable slack force                got=0.000000e+00 want=0.000000e+00
    cable active=0 (want 0)
[PASS] pretension force                 got=1.000000e+05 want=1.000000e+05
[PASS] transverse disp (geom)           got=1.000000e-03 want=1.000000e-03

ALL PASS (0 failures)
[stderr]

Shell cwd was reset to /opt/truss2d
All analytic cases pass, including the geometric transverse-stiffness case (exact `P·L/T0`). Let me do an independent equilibrium-residual check on a richer model (strut, skew springs, geometric coupling) to confirm the converged state truly satisfies equilibrium including the initial-stress transverse forces and reactions.

Equilibrium residual check including geometric forces

shell
cat > /tmp/equil.cpp <<'EOF'
#include <cstdio>
#include <cmath>
#include <vector>
#include "truss2d/model.hpp"
#include "truss2d/solver.hpp"
using namespace truss2d;
int fails=0;

// Check full equilibrium: for each free DOF, sum of internal axial member
// forces + transverse geometric forces + spring reactions + applied loads = 0.
void check_equilibrium(const char* name, const Model& m){
  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;
  // spring reactions (already computed as force on structure)
  for(std::size_t i=0;i<n;i++) resid[i]+=r.reactions[i];
  // member forces: axial + geometric transverse, using converged u and N
  for(std::size_t e=0;e<m.num_elements();e++){
    if(!r.active[e]) continue;
    const Element&el=m.elements()[e];
    Vec2 d=m.element_direction(e); double L=m.element_length(e);
    double c=d.x,s=d.y;
    std::size_t dof[4]={2*el.n1,2*el.n1+1,2*el.n2,2*el.n2+1};
    double N=r.axial_forces[e];
    // axial internal force vector that the member exerts on nodes = -B^T * N
    // (equilibrium: external + (-B^T N) + ... ; here we add member force ON nodes)
    double B[4]={-c,-s,c,s};
    // transverse: relative normal displacement
    double Bn[4]={s,-c,-s,c};
    double un=0; for(int i=0;i<4;i++) un+=Bn[i]*r.displacements[dof[i]];
    double ft=(N/L)*un; // transverse "force" magnitude along Bn
    for(int i=0;i<4;i++){
      // member exerts -B^T N - Bn^T*(N/L)*un on the nodes
      resid[dof[i]] += -B[i]*N - Bn[i]*ft;
    }
  }
  double mx=0; for(double v:resid) mx=std::max(mx,std::fabs(v));
  // scale by load magnitude
  double scale=1.0; for(auto&ld:m.loads()) scale=std::max(scale,std::fabs(ld.value));
  bool ok=mx<=1e-6*scale;
  printf("[%s] %-28s max residual=%.3e (scale %.1e)\n",ok?"PASS":"FAIL",name,mx,scale);
  if(!ok) fails++;
}

int main(){
  double E=2e11,A=1e-3,big=1e15;
  // Strut: compression-only prop. node1 pushed into node0 -> 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,big);m.add_spring(0,0,1,big);m.add_spring(1,0,1,big);
    m.add_load(2,-2e4); // push toward node0 -> compression -> strut active
    auto r=StaticSolver(m).solve();
    printf("    strut force=%.4e active=%d (want compression, active)\n",r.axial_forces[0],(int)r.active[0]);
    if(r.axial_forces[0]>=0||!r.active[0]) fails++;
    check_equilibrium("strut-compress",m);
  }
  // Strut in tension condition -> separates (inactive)
  {
    Model m; m.add_node(0,0); m.add_node(2,0);
    m.add_element(0,1,Section{E,A},MemberKind::Strut);
    double ksm=1e6;
    m.add_spring(0,1,0,big);m.add_spring(0,0,1,big);m.add_spring(1,0,1,big);
    m.add_spring(1,1,0,ksm);
    m.add_load(2,2e4); // pull away -> tension -> strut separates
    auto r=StaticSolver(m).solve();
    printf("    strut(sep) force=%.4e active=%d (want 0, inactive)\n",r.axial_forces[0],(int)r.active[0]);
    if(r.axial_forces[0]!=0.0||r.active[0]) fails++;
  }
  // Skew springs: two nonparallel springs at a node (2D elastic support) + inclined member.
  {
    Model m; m.add_node(0,0); m.add_node(3,4);
    m.add_element(0,1,Section{E,A});
    m.add_spring(0,1,0,big);m.add_spring(0,0,1,big);
    m.add_spring(1,1,0,2e8); m.add_spring(1,1,1,3e8); // skew, nonparallel
    m.add_load(2,5e4); m.add_load(3,-3e4);
    check_equilibrium("skew-springs-incline",m);
  }
  // Single skew roller (one spring) with settlement, inclined member.
  {
    Model m; m.add_node(0,0); m.add_node(3,4);
    m.add_element(0,1,Section{E,A},MemberKind::Bar,2e-4,1e-5,30); // prestrain+thermal
    m.add_spring(0,1,0,big);m.add_spring(0,0,1,big);
    m.add_spring(1,0.6,0.8,1e8,5e-4); // skew roller along member dir, settling
    m.add_spring(1,-0.8,0.6,4e8); // perpendicular finite support
    m.add_load(3,-2e4);
    check_equilibrium("skew-roller-settle-prestrain",m);
  }
  // two_bar example file analog (cable slack + geometric)
  {
    Model m;
    m.add_node(-3,4);m.add_node(3,4);m.add_node(0,0);m.add_node(0,-4);
    m.add_element(2,0,Section{E,A},MemberKind::Cable);
    m.add_element(2,1,Section{E,A},MemberKind::Cable);
    m.add_element(2,3,Section{E,A},MemberKind::Bar);
    double k=5e10;
    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);
    m.add_spring(3,1,0,k);m.add_spring(3,0,1,k);
    m.add_load(4,4e4);
    check_equilibrium("two_bar-example",m);
  }
  printf("\n%s (%d failures)\n", fails?"FAILURES":"ALL PASS", fails);
  return fails?1:0;
}
EOF
g++ -std=c++17 -I/opt/truss2d/include /tmp/equil.cpp /opt/truss2d/build/libtruss2d.a -o /tmp/equil && /tmp/equil
output
strut force=-2.0000e+04 active=1 (want compression, active)
[PASS] strut-compress               max residual=3.638e-12 (scale 2.0e+04)
    strut(sep) force=0.0000e+00 active=0 (want 0, inactive)
[PASS] skew-springs-incline         max residual=1.364e-11 (scale 5.0e+04)
[PASS] skew-roller-settle-prestrain max residual=5.821e-11 (scale 2.0e+04)
[PASS] two_bar-example              max residual=2.910e-11 (scale 4.0e+04)

ALL PASS (0 failures)

[stdout]
    strut force=-2.0000e+04 active=1 (want compression, active)
[PASS] strut-compress               max residual=3.638e-12 (scale 2.0e+04)
    strut(sep) force=0.0000e+00 active=0 (want 0, inactive)
[PASS] skew-springs-incline         max residual=1.364e-11 (scale 5.0e+04)
[PASS] skew-roller-settle-prestrain max residual=5.821e-11 (scale 2.0e+04)
[PASS] two_bar-example              max residual=2.910e-11 (scale 4.0e+04)

ALL PASS (0 failures)
All equilibrium residuals are at machine precision across struts, skew supports, settlement, prestrain+thermal, and unilateral slackening. Let me do final checks on the helper outputs directly, `validate()` propagation, and utilization sign.

Verify helper functions, utilization, validate

shell
cat > /tmp/helpers.cpp <<'EOF'
#include <cstdio>
#include <cmath>
#include "truss2d/model.hpp"
#include "truss2d/solver.hpp"
using namespace truss2d;
int fails=0;
void ck(const char*n,bool c){printf("[%s] %s\n",c?"PASS":"FAIL",n);if(!c)fails++;}
int main(){
  double E=2e11,A=1e-3,big=1e15;
  // helper: element_stiffness_global for horizontal bar
  {
    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,big);m.add_spring(0,0,1,big);m.add_spring(1,0,1,big);
    auto K=StaticSolver(m).element_stiffness_global(0);
    double k=E*A/2.0;
    ck("Ke[0,0]=k",std::fabs(K(0,0)-k)<1e-3);
    ck("Ke[0,2]=-k",std::fabs(K(0,2)+k)<1e-3);
    ck("Ke[1,1]=0 (horizontal)",std::fabs(K(1,1))<1e-6);
    // assemble symmetry
    std::vector<bool> act{true};
    auto S=StaticSolver(m).assemble_system(act);
    bool sym=true; for(std::size_t i=0;i<4;i++)for(std::size_t j=0;j<4;j++) if(std::fabs(S(i,j)-S(j,i))>1e-3)sym=false;
    ck("assemble symmetric",sym);
    // assemble has NO geometric term: node1 y stiffness = spring big only
    ck("assemble no geom (node1 y = big)",std::fabs(S(3,3)-big)<1e-2*big);
    // inactive -> no element contribution
    std::vector<bool> off{false};
    auto S0=StaticSolver(m).assemble_system(off);
    ck("inactive elem dropped",std::fabs(S0(0,0)-big)<1e-2*big);
  }
  // prestrain_force sign: too-long pushes ends apart
  {
    Model m; m.add_node(0,0); m.add_node(2,0);
    double pre=1e-4;
    m.add_element(0,1,Section{E,A},MemberKind::Bar,pre);
    m.add_spring(0,1,0,big);m.add_spring(0,0,1,big);m.add_spring(1,0,1,big);
    std::vector<bool> act{true};
    auto f=StaticSolver(m).prestrain_force(act);
    double g=E*A*pre;
    ck("prestrain f[n1.x]=-g (pushed back)",std::fabs(f[0]+g)<1e-3);
    ck("prestrain f[n2.x]=+g (pushed fwd)",std::fabs(f[2]-g)<1e-3);
  }
  // utilization sign
  {
    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,big);m.add_spring(0,0,1,big);m.add_spring(1,0,1,big);
    m.add_load(2,2.5e8*A); // N=capacity -> util=1
    auto r=StaticSolver(m).solve();
    ck("utilization=+1 at capacity",std::fabs(r.utilization[0]-1.0)<1e-6);
  }
  // validate propagation: no springs -> throws
  {
    Model m; m.add_node(0,0); m.add_node(2,0); m.add_element(0,1,Section{E,A});
    bool threw=false; try{StaticSolver(m).solve();}catch(const std::exception&){threw=true;}
    ck("solve() calls validate (no springs throws)",threw);
  }
  printf("\n%s (%d failures)\n",fails?"FAILURES":"ALL PASS",fails);
  return fails?1:0;
}
EOF
g++ -std=c++17 -I/opt/truss2d/include /tmp/helpers.cpp /opt/truss2d/build/libtruss2d.a -o /tmp/helpers && /tmp/helpers
output
[PASS] Ke[0,0]=k
[PASS] Ke[0,2]=-k
[PASS] Ke[1,1]=0 (horizontal)
[PASS] assemble symmetric
[PASS] assemble no geom (node1 y = big)
[PASS] inactive elem dropped
[PASS] prestrain f[n1.x]=-g (pushed back)
[PASS] prestrain f[n2.x]=+g (pushed fwd)
[PASS] utilization=+1 at capacity
[PASS] solve() calls validate (no springs throws)

ALL PASS (0 failures)

[stdout]
[PASS] Ke[0,0]=k
[PASS] Ke[0,2]=-k
[PASS] Ke[1,1]=0 (horizontal)
[PASS] assemble symmetric
[PASS] assemble no geom (node1 y = big)
[PASS] inactive elem dropped
[PASS] prestrain f[n1.x]=-g (pushed back)
[PASS] prestrain f[n2.x]=+g (pushed fwd)
[PASS] utilization=+1 at capacity
[PASS] solve() calls validate (no springs throws)

ALL PASS (0 failures)
Everything passes. Let me run the official build + ctest one final time to confirm the shipped suite is green, and clean up temp files.

Final ctest run and cleanup

shell
cd /opt/truss2d/build && ctest --output-on-failure 2>&1 | tail -6; rm -f /tmp/verify.cpp /tmp/verify /tmp/equil.cpp /tmp/equil /tmp/helpers.cpp /tmp/helpers
output
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]
    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
I've implemented the four solver functions in `/opt/truss2d/src/solver.cpp`. The build is clean and all tests pass. ## What I implemented **`element_stiffness_global(e)`** , the ordinary elastic 4×4 axial stiffness `(EA/L)·BᵀB` in the global frame, with `B = [-c, -s, c, s]`. Axial only; no kind/prestrain/geometry. **`assemble_system(active)`** , scatters each active member's elastic stiffness into the global matrix and adds each spring's `κ·d̂d̂ᵀ` block at its node. Ordinary elastic + spring only (no initial-stress term, per the contract). **`prestrain_force(active)`** , for each active member, combines mechanical + thermal strain to first order `ε₀ = prestrain + α·dT`, and applies the installed nodal force `EA·ε₀` along `Bᵀ` (a too-long member pushes its ends apart → compressive at zero displacement). **`solve()`** , the full nonlinear analysis: - Calls `Model::validate()` first (exceptions propagate). - Builds the load + settlement RHS (settlement adds `κ·settlement·d̂`; reaction uses stretch relative to that offset). - **Geometric/initial-stress** tangent `K_g = (N/L)·BₙᵀBₙ` with normal `(-s, c)` is added only inside the equilibrium solve , tension stiffens, compression softens, vanishing for rigid translation and axial motion, and never altering the recovered axial force. - A **fixed-point** loop drives the assumed member forces to coincide with the forces recovered from the displacements they produce. The assumed force is seeded with the installed force `-EA·ε₀`, so a transversely-loaded pretensioned member has transverse stiffness from the first solve (a zero seed would make it singular). - An **active-set** loop drops slack cables / separated struts based on each member's trial total axial force, reseeding any reactivated member. - Populates displacements, reactions (spring force on the structure), tension-positive axial forces (exactly 0 for dropped members), and signed utilization `N/(2.5e8·A)`. ## Verification Beyond the shipped smoke test, I checked against independent analytic references (all at machine precision, then removed): axial response, thermal/prestrain sign conventions, spring settlement (displacement, force, and reaction), cable slackening and strut separation/engagement, and the geometric transverse stiffness of a pretensioned member (`v = P·L/T₀` exactly). I also confirmed full equilibrium residuals , including the initial-stress transverse forces and reactions , vanish to ~1e-11 across strut, skew-spring, settlement+prestrain, and the slackening two-bar example; and that the helper outputs (symmetry, no geometric term, inactive-member dropping, prestrain sign) and `validate()` propagation behave correctly. Only `src/solver.cpp` was modified; no headers or signatures changed.
Final diff , what the agent changed

No file changes, the agent left the workspace untouched.

Verifier output , the proof behind the reward
tests/test.sh
== 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.KHSgDv8rtg/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_unilateral.dir/test_unilateral.cpp.o
[ 42%] Building CXX object tests/CMakeFiles/test_core.dir/test_core.cpp.o
[ 52%] Building CXX object tests/CMakeFiles/test_settlement.dir/test_settlement.cpp.o
[ 57%] Building CXX object tests/CMakeFiles/test_geometric.dir/test_geometric.cpp.o
[ 47%] Building CXX object tests/CMakeFiles/test_soak.dir/test_soak.cpp.o
[ 63%] Building CXX object tests/CMakeFiles/test_degenerate.dir/test_degenerate.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_soak
[ 84%] Linking CXX executable test_settlement
[ 84%] Built target test_soak
[ 84%] Built target test_settlement
[ 89%] Linking CXX executable test_unilateral
[ 94%] Linking CXX executable test_core
[ 94%] Built target test_core
[ 94%] Built target test_unilateral
[100%] Linking CXX executable test_geometric
[100%] Built target test_geometric
== Running hidden tests ==
Test project /tmp/tmp.KHSgDv8rtg/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 ...................   Passed    0.00 sec
    Start 4: test_soak
4/6 Test #4: test_soak ........................   Passed    0.00 sec
    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

100% tests passed, 0 tests failed out of 6

Label Time Summary:
hidden    =   0.02 sec*proc (6 tests)

Total Test time (real) =   0.02 sec
PASS: all hidden tests passed

Reproduce this trial: git checkout 2f94510 && PYTHONPATH=src python3 scripts/build_site.py , then open trial/trial_6d9b7f1bf0e04651. Re-running the agent live requires EVAL_PLATFORM_ENABLE_OAUTH_SMOKE=1 and is non-deterministic.

Trial trial_6d9b7f1bf0e04651 · verifier authoritative; classifier explanatory.