truss2d-solver

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

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

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

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

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

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

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:
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
---INSTRUCTION---

[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
---INSTRUCTION---

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

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

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

[ 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/solver.cpp.o
[ 44%] Building CXX object CMakeFiles/truss2d.dir/src/model.cpp.o
[ 55%] Linking CXX static library libtruss2d.a
[ 55%] Built target truss2d
[ 66%] Building CXX object CMakeFiles/truss2d_cli.dir/src/main.cpp.o
[ 77%] Building CXX object tests/CMakeFiles/test_smoke.dir/test_smoke.cpp.o
[ 88%] Linking CXX executable truss2d_cli
[ 88%] Built target truss2d_cli
[100%] Linking CXX executable test_smoke
[100%] Built target test_smoke
Test project /opt/truss2d/build
    Start 1: smoke
1/1 Test #1: smoke ............................   Passed    0.00 sec

100% tests passed, 0 tests failed out of 1

Total Test time (real) =   0.00 sec

[stdout]
[ 11%] Building CXX object CMakeFiles/truss2d.dir/src/dense_matrix.cpp.o
[ 22%] Building CXX object CMakeFiles/truss2d.dir/src/io.cpp.o
[ 33%] Building CXX object CMakeFiles/truss2d.dir/src/solver.cpp.o
[ 44%] Building CXX object CMakeFiles/truss2d.dir/src/model.cpp.o
[ 55%] Linking CXX static library libtruss2d.a
[ 55%] Built target truss2d
[ 66%] Building CXX object CMakeFiles/truss2d_cli.dir/src/main.cpp.o
[ 77%] Building CXX object tests/CMakeFiles/test_smoke.dir/test_smoke.cpp.o
[ 88%] Linking CXX executable truss2d_cli
[ 88%] Built target truss2d_cli
[100%] Linking CXX executable test_smoke
[100%] Built target test_smoke
Test project /opt/truss2d/build
    Start 1: smoke
1/1 Test #1: smoke ............................   Passed    0.00 sec

100% tests passed, 0 tests failed out of 1

Total Test time (real) =   0.00 sec

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 203:  1725 Aborted                 /tmp/verify

[error] tool reported failure

=== case 1 ===
=== case 2 ===
=== case 3 ===
=== case 4 ===
=== case 5 ===
=== case 6 ===
case1 ok N=100000
Shell cwd was reset to /opt/truss2d

[stdout]
=== case 1 ===
=== case 2 ===
=== case 3 ===
=== case 4 ===
=== case 5 ===
=== case 6 ===
case1 ok N=100000
[stderr]

Shell cwd was reset to /opt/truss2d

slack cable N (=0)                         got= 0.000000e+00 want= 0.000000e+00 OK
  cable active=0 (want 0)
slack node1.x u (=F/kb)                    got=-5.000000e-04 want=-5.000000e-04 OK
separated strut N (=0)                     got= 0.000000e+00 want= 0.000000e+00 OK
  strut active=0 (want 0)
separated node1.x u (=F/kb)                got= 5.000000e-04 want= 5.000000e-04 OK
taut cable N (=F)                          got= 4.000000e+04 want= 4.000000e+04 OK
  cable active=1 (want 1)
util (=N/(2.5e8*A))                        got= 1.600000e-01 want= 1.600000e-01 OK

ALL PASS

[stdout]
slack cable N (=0)                         got= 0.000000e+00 want= 0.000000e+00 OK
  cable active=0 (want 0)
slack node1.x u (=F/kb)                    got=-5.000000e-04 want=-5.000000e-04 OK
separated strut N (=0)                     got= 0.000000e+00 want= 0.000000e+00 OK
  strut active=0 (want 0)
separated node1.x u (=F/kb)                got= 5.000000e-04 want= 5.000000e-04 OK
taut cable N (=F)                          got= 4.000000e+04 want= 4.000000e+04 OK
  cable active=1 (want 1)
util (=N/(2.5e8*A))                        got= 1.600000e-01 want= 1.600000e-01 OK

ALL PASS

prestrain N (=-EA*eps0)                      got= 1.000000e+05 want= 1.000000e+05 OK
thermal N (=-EA*al*dT)                       got=-1.920000e+05 want=-1.920000e+05 OK
axial N (=F)                                 got= 5.000000e+04 want= 5.000000e+04 OK
reaction0x (=-F)                             got=-5.000000e+04 want=-5.000000e+04 OK
u1x (=FL/EA)                                 got= 5.000000e-04 want= 5.000000e-04 OK
settlement u1x                               got= 9.900990e-06 want= 9.900990e-06 OK
settlement reaction1x                        got= 9.900990e+02 want= 9.900990e+02 OK
settlement bar N (=k*u)                      got= 9.900990e+02 want= 9.900990e+02 OK
taut-string transverse v                     got= 5.000000e-02 want= 5.000000e-02 OK
taut-string axial N (=T0)                    got= 2.000000e+04 want= 2.000000e+04 OK
string seg2 N == seg1                        got= 2.000000e+04 want= 2.000000e+04 OK

ALL PASS

[stdout]
prestrain N (=-EA*eps0)                      got= 1.000000e+05 want= 1.000000e+05 OK
thermal N (=-EA*al*dT)                       got=-1.920000e+05 want=-1.920000e+05 OK
axial N (=F)                                 got= 5.000000e+04 want= 5.000000e+04 OK
reaction0x (=-F)                             got=-5.000000e+04 want=-5.000000e+04 OK
u1x (=FL/EA)                                 got= 5.000000e-04 want= 5.000000e-04 OK
settlement u1x                               got= 9.900990e-06 want= 9.900990e-06 OK
settlement reaction1x                        got= 9.900990e+02 want= 9.900990e+02 OK
settlement bar N (=k*u)                      got= 9.900990e+02 want= 9.900990e+02 OK
taut-string transverse v                     got= 5.000000e-02 want= 5.000000e-02 OK
taut-string axial N (=T0)                    got= 2.000000e+04 want= 2.000000e+04 OK
string seg2 N == seg1                        got= 2.000000e+04 want= 2.000000e+04 OK

ALL PASS

2-bar N AC                               got= 7.068570e+04 want= 7.071068e+04 OK
2-bar N BC                               got= 7.068570e+04 want= 7.071068e+04 OK
reaction A y (=P/2)                      got= 5.000000e+04 want= 5.000000e+04 OK
reaction B y (=P/2)                      got= 5.000000e+04 want= 5.000000e+04 OK
global sum Fx                            got= 0.000000e+00 want= 0.000000e+00 OK
global sum Fy                            got=-2.910383e-11 want= 0.000000e+00 OK

ALL PASS

[stdout]
2-bar N AC                               got= 7.068570e+04 want= 7.071068e+04 OK
2-bar N BC                               got= 7.068570e+04 want= 7.071068e+04 OK
reaction A y (=P/2)                      got= 5.000000e+04 want= 5.000000e+04 OK
reaction B y (=P/2)                      got= 5.000000e+04 want= 5.000000e+04 OK
global sum Fx                            got= 0.000000e+00 want= 0.000000e+00 OK
global sum Fy                            got=-2.910383e-11 want= 0.000000e+00 OK

ALL PASS

[ 77%] Built target truss2d_cli
[100%] Built target test_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]
[ 77%] Built target truss2d_cli
[100%] Built target test_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