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_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	#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

src/linalg.cpp
src/main.cpp
src/beam.cpp
src/solver.cpp
src/io.cpp
tests/test_smoke.cpp

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

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	#include "beamfd/linalg.hpp"
2	
3	#include <cmath>
4	
5	namespace beamfd {
6	
7	std::vector<double> solve_general(const DenseMatrix& A,
8	                                  const std::vector<double>& b) {
9	    const std::size_t n = A.rows();
10	    if (A.cols() != n) {
11	        throw std::runtime_error("solve_general: matrix is not square");
12	    }
13	    if (b.size() != n) {
14	        throw std::runtime_error("solve_general: right-hand side size mismatch");
15	    }
16	
17	    // Work on mutable copies (Gaussian elimination with partial pivoting).
18	    DenseMatrix M = A;
19	    std::vector<double> x = b;
20	
21	    for (std::size_t col = 0; col < n; ++col) {
22	        // Find the pivot row (largest magnitude in this column at/below the
23	        // diagonal) for numerical stability.
24	        std::size_t pivot = col;
25	        double best = std::fabs(M(col, col));
26	        for (std::size_t r = col + 1; r < n; ++r) {
27	            const double v = std::fabs(M(r, col));
28	            if (v > best) {
29	                best = v;
30	                pivot = r;
31	            }
32	        }
33	        if (best <= 1e-300) {
34	            throw std::runtime_error("solve_general: matrix is singular");
35	        }
36	
37	        // Swap pivot row into place (in both M and the RHS).
38	        if (pivot != col) {
39	            for (std::size_t c = 0; c < n; ++c) {
40	                std::swap(M(col, c), M(pivot, c));
41	            }
42	            std::swap(x[col], x[pivot]);
43	        }
44	
45	        // Eliminate entries below the pivot.
46	        const double diag = M(col, col);
47	        for (std::size_t r = col + 1; r < n; ++r) {
48	            const double factor = M(r, col) / diag;
49	            if (factor == 0.0) continue;
50	            M(r, col) = 0.0;
51	            for (std::size_t c = col + 1; c < n; ++c) {
52	                M(r, c) -= factor * M(col, c);
53	            }
54	            x[r] -= factor * x[col];
55	        }
56	    }
57	
58	    // Back substitution.
59	    for (std::size_t ii = 0; ii < n; ++ii) {
60	        const std::size_t i = n - 1 - ii;
61	        double s = x[i];
62	        for (std::size_t c = i + 1; c < n; ++c) {
63	            s -= M(i, c) * x[c];
64	        }
65	        x[i] = s / M(i, i);
66	    }
67	
68	    return x;
69	}
70	
71	} // namespace beamfd
72

CMakeLists.txt
examples
include
src
tests
---tests---
CMakeLists.txt
test_smoke.cpp
test_util.hpp
---data---
cantilever_tip.beam

[stdout]
CMakeLists.txt
examples
include
src
tests
---tests---
CMakeLists.txt
test_smoke.cpp
test_util.hpp
---data---
cantilever_tip.beam

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