beam-deflection-solver

1	#ifndef BEAMFD_SOLVER_HPP
2	#define BEAMFD_SOLVER_HPP
3	
4	#include <cstddef>
5	#include <vector>
6	
7	#include "beamfd/beam.hpp"
8	#include "beamfd/linalg.hpp"
9	
10	namespace beamfd {
11	
12	/// The assembled finite-difference linear system A w = b for a beam.
13	/// `A` has size num_nodes x num_nodes and `b` has length num_nodes; the
14	/// unknown w is the nodal transverse deflection field. Exposed so tests can
15	/// check the discrete equilibrium residual directly.
16	struct BeamSystem {
17	    DenseMatrix A;
18	    std::vector<double> b;
19	};
20	
21	/// Result of a static beam-bending analysis.
22	struct BeamResult {
23	    /// Nodal transverse deflection w_i [m], length == beam.num_nodes(),
24	    /// ordered from x = 0 (i = 0) to x = L (i = num_intervals).
25	    std::vector<double> deflection;
26	};
27	
28	/// Static bending solver for an Euler-Bernoulli beam with (possibly) variable
29	/// flexural rigidity and an optional elastic foundation. The continuum
30	/// boundary-value problem and the support/load semantics are specified in
31	/// instruction.md; choosing and assembling a suitable finite-difference
32	/// discretization of it, and solving the resulting linear system, is the work
33	/// of this class.
34	class BeamSolver {
35	public:
36	    explicit BeamSolver(const Beam& beam) : beam_(beam) {}
37	
38	    /// Run the full analysis. Calls Beam::validate() first and propagates any
39	    /// exception it throws. Throws std::runtime_error if the assembled system
40	    /// is singular.
41	    BeamResult solve() const;
42	
43	    /// Assemble the finite-difference system A w = b. Exposed for testing and
44	    /// reuse (the discrete equilibrium-residual invariant checks A w - b ~ 0 at
45	    /// the interior nodes).
46	    BeamSystem assemble_system() const;
47	
48	    // ---- Recovery helpers (fully implemented; operate on a deflection field).
49	    // These let the caller recompute physical quantities directly from a
50	    // candidate deflection vector, independent of how it was produced.
51	
52	    /// Bending moment M_i = EI(x_i) * w''(x_i) [N*m] via a central second
53	    /// difference at interior nodes and a one-sided second difference at the
54	    /// ends. `w` must have length num_nodes(). Throws std::invalid_argument
55	    /// otherwise.
56	    std::vector<double> bending_moment_field(const std::vector<double>& w) const;
57	
58	    /// Shear force V_i = (EI w'')'(x_i) [N] via finite differences of the
59	    /// nodal bending moment. `w` must have length num_nodes(). Throws
60	    /// std::invalid_argument otherwise.
61	    std::vector<double> shear_field(const std::vector<double>& w) const;
62	
63	private:
64	    const Beam& beam_;
65	};
66	
67	} // namespace beamfd
68	
69	#endif // BEAMFD_SOLVER_HPP
70

1	#include "beamfd/solver.hpp"
2	
3	#include <cmath>
4	#include <cstddef>
5	#include <stdexcept>
6	#include <vector>
7	
8	#include "beamfd/linalg.hpp"
9	#include "beamfd/types.hpp"
10	
11	namespace beamfd {
12	
13	// ============================================================================
14	// TODO(candidate): Implement the beam-bending finite-difference core.
15	//
16	// Two member functions below are stubs that currently return a zero deflection
17	// field and an empty/zero system, so every analysis is wrong and the grader
18	// fails. Replace the stub bodies of
19	//
20	//     BeamSystem BeamSolver::assemble_system() const;
21	//     BeamResult BeamSolver::solve()          const;
22	//
23	// with correct implementations of the boundary-value problem specified in
24	// instruction.md. Do NOT change the public signatures declared in
25	// include/beamfd/solver.hpp, do NOT modify the recovery helpers below, and
26	// reuse the existing support layer (DenseMatrix, solve_general, and the Beam
27	// accessors EI_at / k_at / q_at / end_moment / applied_end_shear / left /
28	// right / dx / num_intervals / num_nodes / validate).
29	//
30	// The physics, the support and load semantics, the contract, and the I/O
31	// format are described in instruction.md. The required behaviour is enforced
32	// by a hidden grading suite (the shipped tests/ are only a smoke check); design
33	// a discretization that is at least second-order accurate and assemble the
34	// linear system A w = b for the nodal deflection, then solve it with
35	// solve_general. Call Beam::validate() at the start of solve() and let its
36	// exceptions (and any singular-system exception from solve_general) propagate.
37	// ============================================================================
38	
39	BeamSystem BeamSolver::assemble_system() const {
40	    // TODO(candidate): build and return the finite-difference system A w = b.
41	    const std::size_t n = beam_.num_nodes();
42	    BeamSystem sys;
43	    sys.A = DenseMatrix(n, n);  // all zeros -> singular / wrong
44	    sys.b.assign(n, 0.0);
45	    return sys;
46	}
47	
48	BeamResult BeamSolver::solve() const {
49	    // TODO(candidate): validate, assemble, solve, return the deflection field.
50	    const std::size_t n = beam_.num_nodes();
51	    BeamResult result;
52	    result.deflection.assign(n, 0.0);  // trivially zero -> wrong
53	    return result;
54	}
55	
56	// ----------------------------------------------------------------------------
57	// Recovery helpers (already implemented; do NOT modify). These differentiate a
58	// given deflection field so callers can recompute physical quantities (bending
59	// moment, shear) directly from a candidate result, honoring a variable EI(x).
60	
61	namespace {
62	
63	// Second derivative w''(x_i) by central difference (interior) / one-sided
64	// 2nd-order difference (ends).
65	double second_derivative(const std::vector<double>& w, std::size_t i, double h2) {
66	    const std::size_t n = w.size();
67	    if (i == 0) {
68	        return (2.0 * w[0] - 5.0 * w[1] + 4.0 * w[2] - w[3]) / h2;
69	    }
70	    if (i == n - 1) {
71	        return (2.0 * w[n - 1] - 5.0 * w[n - 2] + 4.0 * w[n - 3] - w[n - 4]) / h2;
72	    }
73	    return (w[i - 1] - 2.0 * w[i] + w[i + 1]) / h2;
74	}
75	
76	} // namespace
77	
78	std::vector<double> BeamSolver::bending_moment_field(
79	    const std::vector<double>& w) const {
80	    const std::size_t n = beam_.num_nodes();
81	    if (w.size() != n) {
82	        throw std::invalid_argument("bending_moment_field: deflection size mismatch");
83	    }
84	    const double h2 = beam_.dx() * beam_.dx();
85	    std::vector<double> M(n, 0.0);
86	    for (std::size_t i = 0; i < n; ++i) {
87	        M[i] = beam_.EI_at(i) * second_derivative(w, i, h2);
88	    }
89	    return M;
90	}
91	
92	std::vector<double> BeamSolver::shear_field(const std::vector<double>& w) const {
93	    const std::size_t n = beam_.num_nodes();
94	    if (w.size() != n) {
95	        throw std::invalid_argument("shear_field: deflection size mismatch");
96	    }
97	    if (n < 5) {
98	        throw std::invalid_argument("shear_field: need at least 5 nodes");
99	    }
100	    const double h = beam_.dx();
101	    const double h2 = h * h;
102	    std::vector<double> m(n, 0.0);
103	    for (std::size_t j = 0; j < n; ++j) {
104	        m[j] = beam_.EI_at(j) * second_derivative(w, j, h2);
105	    }
106	    std::vector<double> V(n, 0.0);
107	    for (std::size_t i = 0; i < n; ++i) {
108	        if (i == 0) {
109	            V[i] = (-3.0 * m[0] + 4.0 * m[1] - m[2]) / (2.0 * h);
110	        } else if (i == n - 1) {
111	            V[i] = (3.0 * m[n - 1] - 4.0 * m[n - 2] + m[n - 3]) / (2.0 * h);
112	        } else {
113	            V[i] = (m[i + 1] - m[i - 1]) / (2.0 * h);
114	        }
115	    }
116	    return V;
117	}
118	
119	} // namespace beamfd
120

1	#ifndef BEAMFD_BEAM_HPP
2	#define BEAMFD_BEAM_HPP
3	
4	#include <cstddef>
5	#include <utility>
6	#include <vector>
7	
8	#include "beamfd/types.hpp"
9	
10	namespace beamfd {
11	
12	/// A 1D Euler-Bernoulli beam discretized on a uniform grid, with a possibly
13	/// non-uniform flexural rigidity EI(x), an optional Winkler elastic foundation
14	/// k(x), a (possibly spatially varying) transverse load q(x), and optional
15	/// applied force/moment at free ends.
16	///
17	/// This is a plain data container plus light validation. The numerical core
18	/// (assembling and solving the finite-difference system for the deflection
19	/// field) lives in BeamSolver (see solver.hpp).
20	///
21	/// Geometry / grid:
22	///   - The beam occupies x in [0, L], divided into `num_intervals` equal cells,
23	///     giving num_nodes() = num_intervals + 1 grid points at x_i = i * dx,
24	///     dx = L / num_intervals, for i = 0 .. num_intervals.
25	///
26	/// Fields are stored per node (length num_nodes()):
27	///   - EI_at(i)  flexural rigidity at node i  [N*m^2]      (> 0)
28	///   - k_at(i)   Winkler foundation modulus at node i [N/m^2] (>= 0; 0 = none)
29	///   - q_at(i)   distributed transverse load at node i [N/m]
30	/// plus optional applied actions at free ends (force [N] and moment [N*m]).
31	class Beam {
32	public:
33	    /// Construct a beam of length `length` [m] with a uniform flexural rigidity
34	    /// `EI` [N*m^2], discretized into `num_intervals` equal cells, with the
35	    /// given end supports. EI(x) is initialised constant, k(x) = 0, q(x) = 0.
36	    /// Throws std::invalid_argument if length or EI is non-positive or
37	    /// num_intervals < 2.
38	    Beam(double length, double EI, std::size_t num_intervals, Support left,
39	         Support right);
40	
41	    // ---- Field setters -------------------------------------------------------
42	
43	    /// Set the nodal flexural-rigidity field EI(x_i). Size must equal
44	    /// num_nodes(); every value must be > 0. Throws std::invalid_argument.
45	    void set_ei_nodal(const std::vector<double>& ei);
46	
47	    /// Set EI(x) from piecewise-linear control points (x, value), sampled at
48	    /// each node. Points are taken in the given order; x outside the range is
49	    /// clamped to the nearest endpoint value. Every sampled value must be > 0.
50	    void set_ei_profile(const std::vector<std::pair<double, double>>& points);
51	
52	    /// Set the nodal Winkler foundation field k(x_i) >= 0. Size == num_nodes().
53	    void set_foundation_nodal(const std::vector<double>& k);
54	
55	    /// Set k(x) from piecewise-linear control points (x, value).
56	    void set_foundation_profile(
57	        const std::vector<std::pair<double, double>>& points);
58	
59	    /// Set the nodal distributed-load field q(x_i). Size == num_nodes().
60	    void set_q_nodal(const std::vector<double>& q);
61	
62	    /// Set a uniform distributed load q [N/m] over the whole span (overwrites
63	    /// the load field).
64	    void set_distributed_load(double q);
65	
66	    /// Add a piecewise-linear distributed-load segment ramping from q0 at x0 to
67	    /// q1 at x1 [N/m] to the existing load field. Segments are additive.
68	    void add_load_segment(double x0, double x1, double q0, double q1);
69	
70	    /// Apply a transverse force P [N] at a free end (`at_left_end` -> x = 0,
71	    /// else x = L). Throws std::runtime_error if that end is not Free.
72	    void set_end_load(double P, bool at_left_end);
73	
74	    /// Apply a concentrated moment M [N*m] at a free end. Throws
75	    /// std::runtime_error if that end is not Free.
76	    void set_end_moment(double M, bool at_left_end);
77	
78	    // ---- Accessors -----------------------------------------------------------
79	
80	    double length() const { return length_; }
81	    std::size_t num_intervals() const { return num_intervals_; }
82	    std::size_t num_nodes() const { return num_intervals_ + 1; }
83	    double dx() const { return length_ / static_cast<double>(num_intervals_); }
84	    double node_x(std::size_t i) const;
85	
86	    Support left() const { return left_; }
87	    Support right() const { return right_; }
88	
89	    double EI_at(std::size_t i) const;  ///< nodal flexural rigidity [N*m^2]
90	    double k_at(std::size_t i) const;   ///< nodal foundation modulus [N/m^2]
91	    double q_at(std::size_t i) const;   ///< nodal distributed load [N/m]
92	
93	    /// A representative (maximum nodal) flexural rigidity [N*m^2], used as a
94	    /// scale for residual normalisation. Equals the constant EI for a uniform
95	    /// beam.
96	    double EI() const;
97	
98	    /// True if any nodal foundation modulus is strictly positive.
99	    bool has_foundation() const;
100	
101	    /// Applied end shear V at a free end (`left_end` -> x = 0). A transverse
102	    /// force P at a free end contributes V = -P; 0 if no force is applied.
103	    double applied_end_shear(bool left_end) const;
104	
105	    /// Applied end moment M [N*m] at a free end (0 if none).
106	    double end_moment(bool left_end) const;
107	
108	    /// True if any end force has been applied.
109	    bool has_end_load() const;
110	
111	    /// Throws std::runtime_error if the support / foundation combination cannot
112	    /// suppress rigid-body motion (a free end opposite a non-clamped end, or
113	    /// two free ends, is a mechanism UNLESS an elastic foundation is present).
114	    void validate() const;
115	
116	private:
117	    static std::size_t end_index(bool left_end) { return left_end ? 0u : 1u; }
118	
119	    double length_{0.0};
120	    std::size_t num_intervals_{0};
121	    Support left_{Support::Clamped};
122	    Support right_{Support::Free};
123	    std::vector<double> ei_;  // nodal EI, length num_nodes()
124	    std::vector<double> k_;   // nodal foundation, length num_nodes()
125	    std::vector<double> q_;   // nodal load, length num_nodes()
126	    double force_[2]{0.0, 0.0};   // applied transverse force [left, right]
127	    double moment_[2]{0.0, 0.0};  // applied moment [left, right]
128	};
129	
130	} // namespace beamfd
131	
132	#endif // BEAMFD_BEAM_HPP
133

1	#ifndef BEAMFD_TYPES_HPP
2	#define BEAMFD_TYPES_HPP
3	
4	#include <cstddef>
5	
6	namespace beamfd {
7	
8	/// Type of support condition applied at one end of the beam.
9	///
10	///  - Clamped (built-in / fixed): deflection and slope are both zero
11	///    (w = 0, w' = 0).
12	///  - Pinned (simple support):    deflection is zero and the bending moment
13	///    is zero (w = 0, M = 0).
14	///  - Free:                       the bending moment and shear take their
15	///    applied values (M = M_applied, V = V_applied), both zero when no end
16	///    action is applied.
17	///
18	/// See instruction.md for the continuum boundary-value problem these encode.
19	enum class Support {
20	    Clamped,
21	    Pinned,
22	    Free,
23	};
24	
25	} // namespace beamfd
26	
27	#endif // BEAMFD_TYPES_HPP
28

1	#ifndef BEAMFD_LINALG_HPP
2	#define BEAMFD_LINALG_HPP
3	
4	#include <cstddef>
5	#include <stdexcept>
6	#include <vector>
7	
8	namespace beamfd {
9	
10	/// A small, row-major dense matrix of doubles.
11	///
12	/// Part of the engine's linear-algebra support layer; fully implemented. The
13	/// beam solver uses it to hold the finite-difference system matrix A in
14	/// A w = b. Unlike a structural stiffness matrix, A here is generally NON
15	/// symmetric (boundary-condition rows are one-sided difference stencils), so
16	/// the solve uses a general LU factorization rather than a Cholesky/LDL^T.
17	class DenseMatrix {
18	public:
19	    DenseMatrix() = default;
20	
21	    DenseMatrix(std::size_t rows, std::size_t cols)
22	        : rows_(rows), cols_(cols), data_(rows * cols, 0.0) {}
23	
24	    std::size_t rows() const { return rows_; }
25	    std::size_t cols() const { return cols_; }
26	
27	    double& operator()(std::size_t r, std::size_t c) {
28	        return data_[r * cols_ + c];
29	    }
30	    double operator()(std::size_t r, std::size_t c) const {
31	        return data_[r * cols_ + c];
32	    }
33	
34	    void fill(double v) {
35	        for (auto& x : data_) x = v;
36	    }
37	
38	private:
39	    std::size_t rows_{0};
40	    std::size_t cols_{0};
41	    std::vector<double> data_;
42	};
43	
44	/// Solve the general (possibly non-symmetric) linear system A x = b by LU
45	/// factorization with partial (row) pivoting.
46	///
47	/// `A` is the n-by-n system matrix and `b` is the right-hand side of length n.
48	/// Returns the solution vector x of length n.
49	///
50	/// Throws std::runtime_error if A is not square, if the dimensions are
51	/// inconsistent, or if A is detected to be singular (a pivot column is
52	/// numerically zero), which is how an under-/over-determined boundary-value
53	/// problem surfaces.
54	std::vector<double> solve_general(const DenseMatrix& A,
55	                                  const std::vector<double>& b);
56	
57	} // namespace beamfd
58	
59	#endif // BEAMFD_LINALG_HPP
60

1	#include "beamfd/beam.hpp"
2	
3	#include <algorithm>
4	#include <stdexcept>
5	
6	namespace beamfd {
7	
8	namespace {
9	
10	// Piecewise-linear interpolation of control points (x, value) at coordinate x.
11	// Points are used in the given order (assumed nondecreasing in x); x outside
12	// the range is clamped to the nearest endpoint.
13	double interp(const std::vector<std::pair<double, double>>& pts, double x) {
14	    if (pts.empty()) {
15	        throw std::invalid_argument("Beam: empty profile control points");
16	    }
17	    if (x <= pts.front().first) return pts.front().second;
18	    if (x >= pts.back().first) return pts.back().second;
19	    for (std::size_t s = 1; s < pts.size(); ++s) {
20	        const double x0 = pts[s - 1].first, x1 = pts[s].first;
21	        if (x <= x1) {
22	            const double t = (x1 > x0) ? (x - x0) / (x1 - x0) : 0.0;
23	            return pts[s - 1].second + t * (pts[s].second - pts[s - 1].second);
24	        }
25	    }
26	    return pts.back().second;
27	}
28	
29	} // namespace
30	
31	Beam::Beam(double length, double EI, std::size_t num_intervals, Support left,
32	           Support right)
33	    : length_(length),
34	      num_intervals_(num_intervals),
35	      left_(left),
36	      right_(right) {
37	    if (length_ <= 0.0) {
38	        throw std::invalid_argument("Beam: length must be positive");
39	    }
40	    if (EI <= 0.0) {
41	        throw std::invalid_argument("Beam: EI must be positive");
42	    }
43	    if (num_intervals_ < 2) {
44	        throw std::invalid_argument("Beam: need at least 2 intervals");
45	    }
46	    const std::size_t n = num_nodes();
47	    ei_.assign(n, EI);
48	    k_.assign(n, 0.0);
49	    q_.assign(n, 0.0);
50	}
51	
52	void Beam::set_ei_nodal(const std::vector<double>& ei) {
53	    if (ei.size() != num_nodes()) {
54	        throw std::invalid_argument("Beam::set_ei_nodal: size mismatch");
55	    }
56	    for (double v : ei) {
57	        if (v <= 0.0) {
58	            throw std::invalid_argument("Beam::set_ei_nodal: EI must be positive");
59	        }
60	    }
61	    ei_ = ei;
62	}
63	
64	void Beam::set_ei_profile(const std::vector<std::pair<double, double>>& points) {
65	    std::vector<double> ei(num_nodes());
66	    for (std::size_t i = 0; i < num_nodes(); ++i) ei[i] = interp(points, node_x(i));
67	    set_ei_nodal(ei);
68	}
69	
70	void Beam::set_foundation_nodal(const std::vector<double>& k) {
71	    if (k.size() != num_nodes()) {
72	        throw std::invalid_argument("Beam::set_foundation_nodal: size mismatch");
73	    }
74	    for (double v : k) {
75	        if (v < 0.0) {
76	            throw std::invalid_argument(
77	                "Beam::set_foundation_nodal: k must be non-negative");
78	        }
79	    }
80	    k_ = k;
81	}
82	
83	void Beam::set_foundation_profile(
84	    const std::vector<std::pair<double, double>>& points) {
85	    std::vector<double> k(num_nodes());
86	    for (std::size_t i = 0; i < num_nodes(); ++i) k[i] = interp(points, node_x(i));
87	    set_foundation_nodal(k);
88	}
89	
90	void Beam::set_q_nodal(const std::vector<double>& q) {
91	    if (q.size() != num_nodes()) {
92	        throw std::invalid_argument("Beam::set_q_nodal: size mismatch");
93	    }
94	    q_ = q;
95	}
96	
97	void Beam::set_distributed_load(double q) { q_.assign(num_nodes(), q); }
98	
99	void Beam::add_load_segment(double x0, double x1, double q0, double q1) {
100	    if (x1 < x0) {
101	        std::swap(x0, x1);
102	        std::swap(q0, q1);
103	    }
104	    for (std::size_t i = 0; i < num_nodes(); ++i) {
105	        const double x = node_x(i);
106	        if (x < x0 || x > x1) continue;
107	        const double t = (x1 > x0) ? (x - x0) / (x1 - x0) : 0.0;
108	        q_[i] += q0 + t * (q1 - q0);
109	    }
110	}
111	
112	void Beam::set_end_load(double P, bool at_left_end) {
113	    const Support s = at_left_end ? left_ : right_;
114	    if (s != Support::Free) {
115	        throw std::runtime_error(
116	            "Beam::set_end_load: end load may only be applied at a free end");
117	    }
118	    force_[end_index(at_left_end)] = P;
119	}
120	
121	void Beam::set_end_moment(double M, bool at_left_end) {
122	    const Support s = at_left_end ? left_ : right_;
123	    if (s != Support::Free) {
124	        throw std::runtime_error(
125	            "Beam::set_end_moment: end moment may only be applied at a free end");
126	    }
127	    moment_[end_index(at_left_end)] = M;
128	}
129	
130	double Beam::node_x(std::size_t i) const {
131	    if (i > num_intervals_) {
132	        throw std::out_of_range("Beam::node_x: node index out of range");
133	    }
134	    return static_cast<double>(i) * dx();
135	}
136	
137	double Beam::EI_at(std::size_t i) const {
138	    if (i >= ei_.size()) throw std::out_of_range("Beam::EI_at: index out of range");
139	    return ei_[i];
140	}
141	
142	double Beam::k_at(std::size_t i) const {
143	    if (i >= k_.size()) throw std::out_of_range("Beam::k_at: index out of range");
144	    return k_[i];
145	}
146	
147	double Beam::q_at(std::size_t i) const {
148	    if (i >= q_.size()) throw std::out_of_range("Beam::q_at: index out of range");
149	    return q_[i];
150	}
151	
152	double Beam::EI() const {
153	    double m = 0.0;
154	    for (double v : ei_) m = std::max(m, v);
155	    return m;
156	}
157	
158	bool Beam::has_foundation() const {
159	    for (double v : k_) {
160	        if (v > 0.0) return true;
161	    }
162	    return false;
163	}
164	
165	double Beam::applied_end_shear(bool left_end) const {
166	    return -force_[end_index(left_end)];
167	}
168	
169	double Beam::end_moment(bool left_end) const { return moment_[end_index(left_end)]; }
170	
171	bool Beam::has_end_load() const { return force_[0] != 0.0 || force_[1] != 0.0; }
172	
173	void Beam::validate() const {
174	    const bool left_free = (left_ == Support::Free);
175	    const bool right_free = (right_ == Support::Free);
176	
177	    // An elastic foundation supplies distributed restraint that removes the
178	    // rigid-body modes, so any support combination is well posed when present.
179	    if (has_foundation()) return;
180	
181	    if (left_free && right_free) {
182	        throw std::runtime_error("Beam::validate: free-free beam is a mechanism");
183	    }
184	    if (left_free && right_ != Support::Clamped) {
185	        throw std::runtime_error(
186	            "Beam::validate: a free end requires a clamped opposite end");
187	    }
188	    if (right_free && left_ != Support::Clamped) {
189	        throw std::runtime_error(
190	            "Beam::validate: a free end requires a clamped opposite end");
191	    }
192	}
193	
194	} // namespace beamfd
195

tests/test_util.hpp
tests/test_smoke.cpp
tests/CMakeLists.txt

1	// Smoke test for the beamfd project: a minimal, non-leaky sanity check so the
2	// project compiles and a solver gets quick feedback. It asserts only generic
3	// properties (finite field, pinned ends fixed, a non-trivial deflection, and a
4	// small interior discrete-equilibrium residual) and contains NO closed-form
5	// answer. The authoritative grading is performed by a separate hidden suite.
6	
7	#include <algorithm>
8	#include <cmath>
9	
10	#include "beamfd/beam.hpp"
11	#include "beamfd/solver.hpp"
12	#include "test_util.hpp"
13	
14	using namespace beamfd;
15	
16	BEAMFD_TEST("smoke_pinned_udl_sane") {
17	    Beam beam(2.0, 1000.0, 20, Support::Pinned, Support::Pinned);
18	    beam.set_distributed_load(5.0);
19	    BeamSolver solver(beam);
20	    BeamResult r = solver.solve();
21	    const std::size_t n = beam.num_nodes();
22	
23	    beamfd_test::expect(r.deflection.size() == n, "deflection has num_nodes entries");
24	    for (double w : r.deflection) {
25	        beamfd_test::expect(std::isfinite(w), "deflection is finite");
26	    }
27	    beamfd_test::expect_near(r.deflection[0], 0.0, 1e-9, "left pin fixed");
28	    beamfd_test::expect_near(r.deflection[n - 1], 0.0, 1e-9, "right pin fixed");
29	    beamfd_test::expect(r.deflection[n / 2] > 0.0, "mid-span deflects under load");
30	
31	    // Interior discrete-equilibrium residual A w - b ~ 0 (relative to scale).
32	    BeamSystem sys = solver.assemble_system();
33	    double scale = 0.0;
34	    for (double w : r.deflection) scale = std::max(scale, std::fabs(w));
35	    const double c4 = beam.EI() / std::pow(beam.dx(), 4);
36	    for (std::size_t i = 2; i + 2 < n; ++i) {
37	        double res = -sys.b[i];
38	        for (std::size_t j = 0; j < n; ++j) res += sys.A(i, j) * r.deflection[j];
39	        beamfd_test::expect(std::fabs(res) <= 1e-6 * (c4 * scale) + 1e-9,
40	                            "interior equilibrium residual is small");
41	    }
42	}
43	
44	BEAMFD_TEST_MAIN()
45

1	#ifndef BEAMFD_TEST_UTIL_HPP
2	#define BEAMFD_TEST_UTIL_HPP
3	
4	// Minimal dependency-free test harness for the beamfd suite.
5	//
6	// Each test executable defines test functions and registers them with
7	// BEAMFD_TEST(name) { ... }. main() (provided by BEAMFD_TEST_MAIN) runs every
8	// registered test and returns non-zero if any assertion fails, which CTest
9	// 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 beamfd_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	// Relative-error check (with an absolute floor) for quantities whose scale
56	// varies widely across cases.
57	inline void expect_rel(double got, double want, double rel_tol,
58	                       const std::string& what) {
59	    const double denom = std::fabs(want) > 1e-300 ? std::fabs(want) : 1.0;
60	    if (std::fabs(got - want) / denom > rel_tol) {
61	        throw AssertionError{what + " (got " + std::to_string(got) +
62	                             ", want " + std::to_string(want) +
63	                             ", rel_tol " + std::to_string(rel_tol) + ")"};
64	    }
65	}
66	
67	inline int run_all() {
68	    int failures = 0;
69	    for (const auto& c : registry()) {
70	        try {
71	            c.fn();
72	            std::cout << "[ PASS ] " << c.name << "\n";
73	        } catch (const AssertionError& e) {
74	            std::cout << "[ FAIL ] " << c.name << ": " << e.message << "\n";
75	            ++failures;
76	        } catch (const std::exception& e) {
77	            std::cout << "[ FAIL ] " << c.name
78	                      << ": unexpected exception: " << e.what() << "\n";
79	            ++failures;
80	        } catch (...) {
81	            std::cout << "[ FAIL ] " << c.name << ": unknown exception\n";
82	            ++failures;
83	        }
84	    }
85	    std::cout << "----\n"
86	              << (registry().size() - failures) << "/" << registry().size()
87	              << " tests passed\n";
88	    return failures == 0 ? 0 : 1;
89	}
90	
91	} // namespace beamfd_test
92	
93	#define BEAMFD_CONCAT_INNER(a, b) a##b
94	#define BEAMFD_CONCAT(a, b) BEAMFD_CONCAT_INNER(a, b)
95	
96	#define BEAMFD_TEST(NAME)                                                      \
97	    static void BEAMFD_CONCAT(beamfd_test_fn_, __LINE__)();                    \
98	    static ::beamfd_test::Registrar BEAMFD_CONCAT(beamfd_test_reg_,           \
99	                                                  __LINE__){                   \
100	        NAME, &BEAMFD_CONCAT(beamfd_test_fn_, __LINE__)};                      \
101	    static void BEAMFD_CONCAT(beamfd_test_fn_, __LINE__)()
102	
103	#define BEAMFD_TEST_MAIN()                                                     \
104	    int main() { return ::beamfd_test::run_all(); }
105	
106	#endif // BEAMFD_TEST_UTIL_HPP
107

1	# The shipped project ships only a tiny, non-leaky smoke test using the
2	# header-only harness in test_util.hpp. The authoritative grading is a separate
3	# HIDDEN suite injected by the verifier (tests/test.sh) at grade time, so do
4	# not rely on this smoke test to characterize correctness.
5	
6	add_executable(test_smoke test_smoke.cpp)
7	target_link_libraries(test_smoke PRIVATE beamfd)
8	target_include_directories(test_smoke PRIVATE ${CMAKE_CURRENT_SOURCE_DIR})
9	add_test(NAME smoke COMMAND test_smoke)
10

/opt/beamfd:
CMakeLists.txt
examples
include
src
tests

/opt/beamfd/examples:
cantilever_tip.beam

/opt/beamfd/include:
beamfd

/opt/beamfd/include/beamfd:
beam.hpp
io.hpp
linalg.hpp
solver.hpp
types.hpp

/opt/beamfd/src:
beam.cpp
io.cpp
linalg.cpp
main.cpp
solver.cpp

/opt/beamfd/tests:
CMakeLists.txt
test_smoke.cpp
test_util.hpp

[stdout]
/opt/beamfd:
CMakeLists.txt
examples
include
src
tests

/opt/beamfd/examples:
cantilever_tip.beam

/opt/beamfd/include:
beamfd

/opt/beamfd/include/beamfd:
beam.hpp
io.hpp
linalg.hpp
solver.hpp
types.hpp

/opt/beamfd/src:
beam.cpp
io.cpp
linalg.cpp
main.cpp
solver.cpp

/opt/beamfd/tests:
CMakeLists.txt
test_smoke.cpp
test_util.hpp

BUILD_OK
Test project /opt/beamfd/build
    Start 1: smoke
1/1 Test #1: smoke ............................   Passed    0.00 sec

100% tests passed, 0 tests failed out of 1

Total Test time (real) =   0.00 sec

[stdout]
BUILD_OK
Test project /opt/beamfd/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

clamped-clamped const EI          N= 20 err=3.70e-02  N= 40 err=9.62e-03 (ord 1.94)  N= 80 err=2.45e-03 (ord 1.97)  N=160 err=6.19e-04 (ord 1.99)  N=320 err=1.56e-04 (ord 1.99)
pinned-pinned const EI (sin)      N= 20 err=6.66e-04  N= 40 err=4.32e-04 (ord 0.62)  N= 80 err=1.82e-04 (ord 1.24)  N=160 err=5.49e-05 (ord 1.73)  N=320 err=1.49e-05 (ord 1.88)
clamped-free const EI (x^4)       N= 20 err=2.83e-02  N= 40 err=6.97e-03 (ord 2.02)  N= 80 err=1.73e-03 (ord 2.01)  N=160 err=4.31e-04 (ord 2.00)  N=320 err=1.08e-04 (ord 2.00)
clamped-clamped linear EI         N= 20 err=3.70e-02  N= 40 err=9.63e-03 (ord 1.94)  N= 80 err=2.45e-03 (ord 1.97)  N=160 err=6.19e-04 (ord 1.99)  N=320 err=1.56e-04 (ord 1.99)
pinned-pinned + foundation        N= 20 err=4.14e-04  N= 40 err=2.45e-04 (ord 0.76)  N= 80 err=1.01e-04 (ord 1.28)  N=160 err=3.03e-05 (ord 1.74)  N=320 err=8.20e-06 (ord 1.89)
free-free + foundation (cos)      N= 20 err=2.21e-02  N= 40 err=8.25e-03 (ord 1.42)  N= 80 err=2.41e-03 (ord 1.78)  N=160 err=6.47e-04 (ord 1.90)  N=320 err=1.67e-04 (ord 1.95)

[stdout]
clamped-clamped const EI          N= 20 err=3.70e-02  N= 40 err=9.62e-03 (ord 1.94)  N= 80 err=2.45e-03 (ord 1.97)  N=160 err=6.19e-04 (ord 1.99)  N=320 err=1.56e-04 (ord 1.99)
pinned-pinned const EI (sin)      N= 20 err=6.66e-04  N= 40 err=4.32e-04 (ord 0.62)  N= 80 err=1.82e-04 (ord 1.24)  N=160 err=5.49e-05 (ord 1.73)  N=320 err=1.49e-05 (ord 1.88)
clamped-free const EI (x^4)       N= 20 err=2.83e-02  N= 40 err=6.97e-03 (ord 2.02)  N= 80 err=1.73e-03 (ord 2.01)  N=160 err=4.31e-04 (ord 2.00)  N=320 err=1.08e-04 (ord 2.00)
clamped-clamped linear EI         N= 20 err=3.70e-02  N= 40 err=9.63e-03 (ord 1.94)  N= 80 err=2.45e-03 (ord 1.97)  N=160 err=6.19e-04 (ord 1.99)  N=320 err=1.56e-04 (ord 1.99)
pinned-pinned + foundation        N= 20 err=4.14e-04  N= 40 err=2.45e-04 (ord 0.76)  N= 80 err=1.01e-04 (ord 1.28)  N=160 err=3.03e-05 (ord 1.74)  N=320 err=8.20e-06 (ord 1.89)
free-free + foundation (cos)      N= 20 err=2.21e-02  N= 40 err=8.25e-03 (ord 1.42)  N= 80 err=2.41e-03 (ord 1.78)  N=160 err=6.47e-04 (ord 1.90)  N=320 err=1.67e-04 (ord 1.95)

clamped-clamped const EI          N= 20 err=3.70e-02  N= 40 err=9.62e-03 (ord 1.94)  N= 80 err=2.45e-03 (ord 1.97)  N=160 err=6.19e-04 (ord 1.99)  N=320 err=1.56e-04 (ord 1.99)
pinned-pinned const EI (sin)      N= 20 err=4.12e-03  N= 40 err=1.03e-03 (ord 2.00)  N= 80 err=2.57e-04 (ord 2.00)  N=160 err=6.43e-05 (ord 2.00)  N=320 err=1.61e-05 (ord 2.00)
clamped-free const EI (x^4)       N= 20 err=2.83e-02  N= 40 err=6.97e-03 (ord 2.02)  N= 80 err=1.73e-03 (ord 2.01)  N=160 err=4.31e-04 (ord 2.00)  N=320 err=1.08e-04 (ord 2.00)
clamped-clamped linear EI         N= 20 err=3.70e-02  N= 40 err=9.63e-03 (ord 1.94)  N= 80 err=2.45e-03 (ord 1.97)  N=160 err=6.19e-04 (ord 1.99)  N=320 err=1.56e-04 (ord 1.99)
pinned-pinned + foundation        N= 20 err=2.26e-03  N= 40 err=5.65e-04 (ord 2.00)  N= 80 err=1.41e-04 (ord 2.00)  N=160 err=3.53e-05 (ord 2.00)  N=320 err=8.82e-06 (ord 2.00)
free-free + foundation (cos)      N= 20 err=2.21e-02  N= 40 err=8.25e-03 (ord 1.42)  N= 80 err=2.41e-03 (ord 1.78)  N=160 err=6.47e-04 (ord 1.90)  N=320 err=1.67e-04 (ord 1.95)

[stdout]
clamped-clamped const EI          N= 20 err=3.70e-02  N= 40 err=9.62e-03 (ord 1.94)  N= 80 err=2.45e-03 (ord 1.97)  N=160 err=6.19e-04 (ord 1.99)  N=320 err=1.56e-04 (ord 1.99)
pinned-pinned const EI (sin)      N= 20 err=4.12e-03  N= 40 err=1.03e-03 (ord 2.00)  N= 80 err=2.57e-04 (ord 2.00)  N=160 err=6.43e-05 (ord 2.00)  N=320 err=1.61e-05 (ord 2.00)
clamped-free const EI (x^4)       N= 20 err=2.83e-02  N= 40 err=6.97e-03 (ord 2.02)  N= 80 err=1.73e-03 (ord 2.01)  N=160 err=4.31e-04 (ord 2.00)  N=320 err=1.08e-04 (ord 2.00)
clamped-clamped linear EI         N= 20 err=3.70e-02  N= 40 err=9.63e-03 (ord 1.94)  N= 80 err=2.45e-03 (ord 1.97)  N=160 err=6.19e-04 (ord 1.99)  N=320 err=1.56e-04 (ord 1.99)
pinned-pinned + foundation        N= 20 err=2.26e-03  N= 40 err=5.65e-04 (ord 2.00)  N= 80 err=1.41e-04 (ord 2.00)  N=160 err=3.53e-05 (ord 2.00)  N=320 err=8.82e-06 (ord 2.00)
free-free + foundation (cos)      N= 20 err=2.21e-02  N= 40 err=8.25e-03 (ord 1.42)  N= 80 err=2.41e-03 (ord 1.78)  N=160 err=6.47e-04 (ord 1.90)  N=320 err=1.67e-04 (ord 1.95)

free-free + foundation (cos)      N= 20 err=2.21e-02  N= 40 err=8.25e-03 (ord 1.42)  N= 80 err=2.41e-03 (ord 1.78)  N=160 err=6.47e-04 (ord 1.90)  N=320 err=1.67e-04 (ord 1.95)
clamped-free transcendental       N= 20 err=5.92e-01  N= 40 err=5.86e-01 (ord 0.02)  N= 80 err=5.84e-01 (ord 0.00)  N=160 err=5.84e-01 (ord 0.00)  N=320 err=5.84e-01 (ord 0.00)

[stdout]
free-free + foundation (cos)      N= 20 err=2.21e-02  N= 40 err=8.25e-03 (ord 1.42)  N= 80 err=2.41e-03 (ord 1.78)  N=160 err=6.47e-04 (ord 1.90)  N=320 err=1.67e-04 (ord 1.95)
clamped-free transcendental       N= 20 err=5.92e-01  N= 40 err=5.86e-01 (ord 0.02)  N= 80 err=5.84e-01 (ord 0.00)  N=160 err=5.84e-01 (ord 0.00)  N=320 err=5.84e-01 (ord 0.00)

clamped-free transcendental       N= 20 err=1.49e-02  N= 40 err=3.65e-03 (ord 2.03)  N= 80 err=9.02e-04 (ord 2.02)  N=160 err=2.24e-04 (ord 2.01)  N=320 err=5.58e-05 (ord 2.01)

[stdout]
clamped-free transcendental       N= 20 err=1.49e-02  N= 40 err=3.65e-03 (ord 2.03)  N= 80 err=9.02e-04 (ord 2.02)  N=160 err=2.24e-04 (ord 2.01)  N=320 err=5.58e-05 (ord 2.01)

clamped-pinned const EI           N= 20 err=1.65e-02  N= 40 err=4.25e-03 (ord 1.96)  N= 80 err=1.08e-03 (ord 1.98)  N=160 err=2.72e-04 (ord 1.99)  N=320 err=6.83e-05 (ord 1.99)
clamped-free variable EI (x^4)    N= 20 err=2.93e-02  N= 40 err=7.43e-03 (ord 1.98)  N= 80 err=1.87e-03 (ord 1.99)  N=160 err=4.69e-04 (ord 2.00)  N=320 err=1.18e-04 (ord 2.00)

[stdout]
clamped-pinned const EI           N= 20 err=1.65e-02  N= 40 err=4.25e-03 (ord 1.96)  N= 80 err=1.08e-03 (ord 1.98)  N=160 err=2.72e-04 (ord 1.99)  N=320 err=6.83e-05 (ord 1.99)
clamped-free variable EI (x^4)    N= 20 err=2.93e-02  N= 40 err=7.43e-03 (ord 1.98)  N= 80 err=1.87e-03 (ord 1.99)  N=160 err=4.69e-04 (ord 2.00)  N=320 err=1.18e-04 (ord 2.00)

39	// ----------------------------------------------------------------------------
40	// Finite-difference discretization of
41	//
42	//     d^2/dx^2 ( EI(x) w''(x) ) + k(x) w(x) = q(x)            on x in [0, L]
43	//
44	// on the uniform grid x_i = i*dx, i = 0..N (n = N+1 nodes), with two boundary
45	// conditions imposed at each end. Everything is second-order accurate in dx.
46	//
47	// Interior nodes (2 <= i <= N-2) use the conservative ("divergence form")
48	// stencil: the nodal bending moment M_j = EI_j * w''_j is formed with a central
49	// second difference, and d^2/dx^2 of that field is taken with a second central
50	// difference of M. This double application keeps the operator self-adjoint and
51	// reduces, for constant EI, to the classic EI*(1,-4,6,-4,1)/dx^4 biharmonic.
52	//
53	// The two near-boundary equations at each end (rows 0,1 on the left and N-1,N
54	// on the right) are replaced by that end's two boundary conditions, discretized
55	// to the same (second) order:
56	//   * essential (w = 0)        -> a unit row;
57	//   * clamped slope (w' = 0)   -> a 2nd-order one-sided first difference;
58	//   * pinned/free moment       -> EI * (2nd-order one-sided second difference);
59	//   * free shear (EI w'')'     -> a 2nd-order one-sided first difference of the
60	//                                 nodal moment field M_j (each M_j a central
61	//                                 second difference), so the sampled truncation
62	//                                 error stays smooth and the order is preserved.
63	// ----------------------------------------------------------------------------
64	
65	BeamSystem BeamSolver::assemble_system() const {
66	    const std::size_t n = beam_.num_nodes();
67	    const std::size_t N = beam_.num_intervals();  // last node index
68	    const double h = beam_.dx();
69	    const double h2 = h * h;
70	    const double inv_h2 = 1.0 / h2;
71	    const double inv_h4 = 1.0 / (h2 * h2);
72	
73	    BeamSystem sys;
74	    sys.A = DenseMatrix(n, n);
75	    sys.b.assign(n, 0.0);
76	    DenseMatrix& A = sys.A;
77	
78	    // Adds the central-second-difference moment M_j = EI_j*(w_{j-1}-2w_j+w_{j+1})
79	    // / h^2, weighted by `scale`, into `row`. Node j must be interior so that
80	    // j-1 and j+1 exist.
81	    auto add_moment = [&](std::size_t row, std::size_t j, double scale) {
82	        const double c = scale * beam_.EI_at(j) * inv_h2;
83	        A(row, j - 1) += c;
84	        A(row, j) += -2.0 * c;
85	        A(row, j + 1) += c;
86	    };
87	
88	    // Adds EI_p * w''(x_p), weighted by `scale`, into `row`, using a one-sided
89	    // second difference marching in direction `dir` (+1 from the left end, -1
90	    // from the right). Second order when four points are available, otherwise a
91	    // first-order three-point fallback for very coarse grids.
92	    auto add_end_moment_term = [&](std::size_t row, std::size_t p, int dir,
93	                                   double scale) {
94	        const double c = scale * beam_.EI_at(p) * inv_h2;
95	        const long lp = static_cast<long>(p);
96	        const long last = static_cast<long>(N);
97	        const long p3 = lp + 3L * dir;
98	        if (p3 >= 0 && p3 <= last) {
99	            A(row, p) += 2.0 * c;
100	            A(row, static_cast<std::size_t>(lp + dir)) += -5.0 * c;
101	            A(row, static_cast<std::size_t>(lp + 2 * dir)) += 4.0 * c;
102	            A(row, static_cast<std::size_t>(p3)) += -1.0 * c;
103	        } else {
104	            A(row, p) += c;
105	            A(row, static_cast<std::size_t>(lp + dir)) += -2.0 * c;
106	            A(row, static_cast<std::size_t>(lp + 2 * dir)) += c;
107	        }
108	    };
109	
110	    // --- Interior conservative rows: (M_{i-1} - 2 M_i + M_{i+1})/h^2 + k w = q.
111	    for (std::size_t i = 2; i + 2 <= N; ++i) {
112	        const double EIm = beam_.EI_at(i - 1);
113	        const double EIc = beam_.EI_at(i);
114	        const double EIp = beam_.EI_at(i + 1);
115	        A(i, i - 2) += EIm * inv_h4;
116	        A(i, i - 1) += (-2.0 * EIm - 2.0 * EIc) * inv_h4;
117	        A(i, i) += (EIm + 4.0 * EIc + EIp) * inv_h4 + beam_.k_at(i);
118	        A(i, i + 1) += (-2.0 * EIc - 2.0 * EIp) * inv_h4;
119	        A(i, i + 2) += EIp * inv_h4;
120	        sys.b[i] = beam_.q_at(i);
121	    }
122	
123	    // --- Left boundary: rows 0 (end node) and 1 (adjacent node).
124	    switch (beam_.left()) {
125	        case Support::Clamped:
126	            // w_0 = 0 ; w'(0) = 0  (forward 2nd-order first difference).
127	            A(0, 0) = 1.0;
128	            A(1, 0) = -3.0 / (2.0 * h);
129	            A(1, 1) = 4.0 / (2.0 * h);
130	            A(1, 2) = -1.0 / (2.0 * h);
131	            break;
132	        case Support::Pinned: {
133	            // w_0 = 0 ; zero moment via the central reflection w_{-1} = -w_1
134	            // (so the discrete w''(0) = 0), with the conservative governing
135	            // equation applied at node 1 (M_0 collapses to 0).
136	            A(0, 0) = 1.0;
137	            const double EI1 = beam_.EI_at(1), EI2 = beam_.EI_at(2);
138	            A(1, 1) += (4.0 * EI1 + EI2) * inv_h4 + beam_.k_at(1);
139	            A(1, 2) += (-2.0 * EI1 - 2.0 * EI2) * inv_h4;
140	            A(1, 3) += EI2 * inv_h4;
141	            sys.b[1] = beam_.q_at(1);
142	            break;
143	        }
144	        case Support::Free:
145	            // EI_0 w''(0) = M_applied ; (EI w'')'(0) = V_applied.
146	            add_end_moment_term(0, 0, +1, 1.0);
147	            sys.b[0] = beam_.end_moment(true);
148	            add_moment(1, 1, -5.0 / (2.0 * h));
149	            add_moment(1, 2, 8.0 / (2.0 * h));
150	            add_moment(1, 3, -3.0 / (2.0 * h));
151	            sys.b[1] = beam_.applied_end_shear(true);
152	            break;
153	    }
154	
155	    // --- Right boundary: rows N (end node) and N-1 (adjacent node).
156	    switch (beam_.right()) {
157	        case Support::Clamped:
158	            // w_N = 0 ; w'(L) = 0  (backward 2nd-order first difference).
159	            A(N, N) = 1.0;
160	            A(N - 1, N) = 3.0 / (2.0 * h);
161	            A(N - 1, N - 1) = -4.0 / (2.0 * h);
162	            A(N - 1, N - 2) = 1.0 / (2.0 * h);
163	            break;
164	        case Support::Pinned: {
165	            // w_N = 0 ; zero moment via the reflection w_{N+1} = -w_{N-1}, with
166	            // the conservative governing equation applied at node N-1.
167	            A(N, N) = 1.0;
168	            const double EIa = beam_.EI_at(N - 1), EIb = beam_.EI_at(N - 2);
169	            A(N - 1, N - 1) += (4.0 * EIa + EIb) * inv_h4 + beam_.k_at(N - 1);
170	            A(N - 1, N - 2) += (-2.0 * EIa - 2.0 * EIb) * inv_h4;
171	            A(N - 1, N - 3) += EIb * inv_h4;
172	            sys.b[N - 1] = beam_.q_at(N - 1);
173	            break;
174	        }
175	        case Support::Free:
176	            // EI_N w''(L) = M_applied ; (EI w'')'(L) = V_applied.
177	            add_end_moment_term(N, N, -1, 1.0);
178	            sys.b[N] = beam_.end_moment(false);
179	            add_moment(N - 1, N - 1, 5.0 / (2.0 * h));
180	            add_moment(N - 1, N - 2, -8.0 / (2.0 * h));
181	            add_moment(N - 1, N - 3, 3.0 / (2.0 * h));
182	            sys.b[N - 1] = beam_.applied_end_shear(false);
183	            break;
184	    }
185	
186	    return sys;
187	}
188	
189	BeamResult BeamSolver::solve() const {
190	    beam_.validate();  // propagate mechanism / configuration errors
191	    const BeamSystem sys = assemble_system();
192	    BeamResult result;
193	    result.deflection = solve_general(sys.A, sys.b);  // propagates if singular
194	    return result;
195	}
196	
197	// ----------------------------------------------------------------------------
198	// Recovery helpers (already implemented; do NOT modify). These differentiate a
199	// given deflection field so callers can recompute physical quantities (bending
200	// moment, shear) directly from a candidate result, honoring a variable EI(x).
201	
202	namespace {
203	
204	// Second derivative w''(x_i) by central difference (interior) / one-sided
205	// 2nd-order difference (ends).
206	double second_derivative(const std::vector<double>& w, std::size_t i, double h2) {
207	    const std::size_t n = w.size();
208	    if (i == 0) {
209	        return (2.0 * w[0] - 5.0 * w[1] + 4.0 * w[2] - w[3]) / h2;
210	    }
211	    if (i == n - 1) {
212	        return (2.0 * w[n - 1] - 5.0 * w[n - 2] + 4.0 * w[n - 3] - w[n - 4]) / h2;
213	    }

OK
    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]
OK
    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

build/tests/test_smoke
build/beamfd_cli
# Cantilever beam with a downward point load at the free tip.
# Clamped (built-in) at x = 0, free at x = L.
length 3.0          # span [m]
ei     2.0e6        # flexural rigidity EI [N*m^2]
cells  200          # number of equal finite-difference intervals

# Left end clamped, right end free.
support clamped free

# Downward point load of 1 kN at the free right tip.
endload right 1000.0

[stdout]
build/tests/test_smoke
build/beamfd_cli
# Cantilever beam with a downward point load at the free tip.
# Clamped (built-in) at x = 0, free at x = L.
length 3.0          # span [m]
ei     2.0e6        # flexural rigidity EI [N*m^2]
cells  200          # number of equal finite-difference intervals

# Left end clamped, right end free.
support clamped free

# Downward point load of 1 kN at the free right tip.
endload right 1000.0

# beam analysis report
length=3.000000e+00 EI(max)=2.000000e+06 nodes=201 dx=1.500000e-02

[fields] (i, x, w[m], M=EI*w''[N*m], V=(EI*w'')'[N])
0 0.000000e+00 -1.532131e-11 3.000000e+03 -1.000000e+03
1 1.500000e-02 1.678912e-07 2.985000e+03 -1.000000e+03
2 3.000000e-02 6.716103e-07 2.970000e+03 -1.000000e+03
3 4.500000e-02 1.509454e-06 2.955000e+03 -1.000000e+03
4 6.000000e-02 2.679736e-06 2.940000e+03 -1.000000e+03
5 7.500000e-02 4.180767e-06 2.925000e+03 -1.000000e+03
6 9.000000e-02 6.010862e-06 2.910000e+03 -1.000000e+03
7 1.050000e-01 8.168331e-06 2.895000e+03 -1.000000e+03
8 1.200000e-01 1.065149e-05 2.880000e+03 -1.000000e+03
9 1.350000e-01 1.345864e-05 2.865000e+03 -1.000000e+03
10 1.500000e-01 1.658811e-05 2.850000e+03 -1.000000e+03
11 1.650000e-01 2.003821e-05 2.835000e+03 -1.000000e+03
12 1.800000e-01 2.380724e-05 2.820000e+03 -1.000000e+03
13 1.950000e-01 2.789352e-05 2.805000e+03 -1.000000e+03
14 2.100000e-01 3.229536e-05 2.790000e+03 -1.000000e+03
15 2.250000e-01 3.701108e-05 2.775000e+03 -1.000000e+03
16 2.400000e-01 4.203899e-05 2.760000e+03 -1.000000e+03
17 2.550000e-01 4.737740e-05 2.745000e+03 -1.000000e+03
18 2.700000e-01 5.302461e-05 2.730000e+03 -1.000000e+03
19 2.850000e-01 5.897896e-05 2.715000e+03 -1.000000e+03
20 3.000000e-01 6.523874e-05 2.700000e+03 -1.000000e+03
21 3.150000e-01 7.180227e-05 2.685000e+03 -1.000000e+03
22 3.300000e-01 7.866787e-05 2.670000e+03 -1.000000e+03
23 3.450000e-01 8.583383e-05 2.655000e+03 -1.000000e+03
24 3.600000e-01 9.329849e-05 2.640000e+03 -1.000000e+03
25 3.750000e-01 1.010601e-04 2.625000e+03 -1.000000e+03
26 3.900000e-01 1.091171e-04 2.610000e+03 -1.000000e+03
27 4.050000e-01 1.174677e-04 2.595000e+03 -1.000000e+03
28 4.200000e-01 1.261102e-04 2.580000e+03 -1.000000e+03
29 4.350000e-01 1.350430e-04 2.565000e+03 -1.000000e+03
30 4.500000e-01 1.442644e-04 2.550000e+03 -1.000000e+03
31 4.650000e-01 1.537726e-04 2.535000e+03 -1.000000e+03
32 4.800000e-01 1.635660e-04 2.520000e+03 -1.000000e+03
33 4.950000e-01 1.736429e-04 2.505000e+03 -1.000000e+03
34 5.100000e-01 1.840016e-04 2.490000e+03 -1.000000e+03
35 5.250000e-01 1.946405e-04 2.475000e+03 -1.000000e+03

[stdout]
# beam analysis report
length=3.000000e+00 EI(max)=2.000000e+06 nodes=201 dx=1.500000e-02

[fields] (i, x, w[m], M=EI*w''[N*m], V=(EI*w'')'[N])
0 0.000000e+00 -1.532131e-11 3.000000e+03 -1.000000e+03
1 1.500000e-02 1.678912e-07 2.985000e+03 -1.000000e+03
2 3.000000e-02 6.716103e-07 2.970000e+03 -1.000000e+03
3 4.500000e-02 1.509454e-06 2.955000e+03 -1.000000e+03
4 6.000000e-02 2.679736e-06 2.940000e+03 -1.000000e+03
5 7.500000e-02 4.180767e-06 2.925000e+03 -1.000000e+03
6 9.000000e-02 6.010862e-06 2.910000e+03 -1.000000e+03
7 1.050000e-01 8.168331e-06 2.895000e+03 -1.000000e+03
8 1.200000e-01 1.065149e-05 2.880000e+03 -1.000000e+03
9 1.350000e-01 1.345864e-05 2.865000e+03 -1.000000e+03
10 1.500000e-01 1.658811e-05 2.850000e+03 -1.000000e+03
11 1.650000e-01 2.003821e-05 2.835000e+03 -1.000000e+03
12 1.800000e-01 2.380724e-05 2.820000e+03 -1.000000e+03
13 1.950000e-01 2.789352e-05 2.805000e+03 -1.000000e+03
14 2.100000e-01 3.229536e-05 2.790000e+03 -1.000000e+03
15 2.250000e-01 3.701108e-05 2.775000e+03 -1.000000e+03
16 2.400000e-01 4.203899e-05 2.760000e+03 -1.000000e+03
17 2.550000e-01 4.737740e-05 2.745000e+03 -1.000000e+03
18 2.700000e-01 5.302461e-05 2.730000e+03 -1.000000e+03
19 2.850000e-01 5.897896e-05 2.715000e+03 -1.000000e+03
20 3.000000e-01 6.523874e-05 2.700000e+03 -1.000000e+03
21 3.150000e-01 7.180227e-05 2.685000e+03 -1.000000e+03
22 3.300000e-01 7.866787e-05 2.670000e+03 -1.000000e+03
23 3.450000e-01 8.583383e-05 2.655000e+03 -1.000000e+03
24 3.600000e-01 9.329849e-05 2.640000e+03 -1.000000e+03
25 3.750000e-01 1.010601e-04 2.625000e+03 -1.000000e+03
26 3.900000e-01 1.091171e-04 2.610000e+03 -1.000000e+03
27 4.050000e-01 1.174677e-04 2.595000e+03 -1.000000e+03
28 4.200000e-01 1.261102e-04 2.580000e+03 -1.000000e+03
29 4.350000e-01 1.350430e-04 2.565000e+03 -1.000000e+03
30 4.500000e-01 1.442644e-04 2.550000e+03 -1.000000e+03
31 4.650000e-01 1.537726e-04 2.535000e+03 -1.000000e+03
32 4.800000e-01 1.635660e-04 2.520000e+03 -1.000000e+03
33 4.950000e-01 1.736429e-04 2.505000e+03 -1.000000e+03
34 5.100000e-01 1.840016e-04 2.490000e+03 -1.000000e+03
35 5.250000e-01 1.946405e-04 2.475000e+03 -1.000000e+03