beam-deflection-solver

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

./CMakeLists.txt
./examples/cantilever_tip.beam
./src/beam.cpp
./src/io.cpp
./src/linalg.cpp
./src/main.cpp
./src/solver.cpp
./tests/CMakeLists.txt
./tests/test_smoke.cpp
---TESTS---
CMakeLists.txt
test_smoke.cpp
test_util.hpp

[stdout]
./CMakeLists.txt
./examples/cantilever_tip.beam
./src/beam.cpp
./src/io.cpp
./src/linalg.cpp
./src/main.cpp
./src/solver.cpp
./tests/CMakeLists.txt
./tests/test_smoke.cpp
---TESTS---
CMakeLists.txt
test_smoke.cpp
test_util.hpp

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	# Cantilever beam with a downward point load at the free tip.
2	# Clamped (built-in) at x = 0, free at x = L.
3	length 3.0          # span [m]
4	ei     2.0e6        # flexural rigidity EI [N*m^2]
5	cells  200          # number of equal finite-difference intervals
6	
7	# Left end clamped, right end free.
8	support clamped free
9	
10	# Downward point load of 1 kN at the free right tip.
11	endload right 1000.0
12

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.01 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.01 sec

1	cmake_minimum_required(VERSION 3.16)
2	project(beamfd LANGUAGES CXX)
3	
4	set(CMAKE_CXX_STANDARD 17)
5	set(CMAKE_CXX_STANDARD_REQUIRED ON)
6	set(CMAKE_CXX_EXTENSIONS OFF)
7	
8	if(NOT CMAKE_BUILD_TYPE)
9	  set(CMAKE_BUILD_TYPE Release)
10	endif()
11	
12	add_compile_options(-Wall -Wextra)
13	
14	# --- Core library -----------------------------------------------------------
15	add_library(beamfd
16	  src/linalg.cpp
17	  src/beam.cpp
18	  src/solver.cpp
19	  src/io.cpp
20	)
21	target_include_directories(beamfd PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}/include)
22	
23	# --- Demo CLI ---------------------------------------------------------------
24	add_executable(beamfd_cli src/main.cpp)
25	target_link_libraries(beamfd_cli PRIVATE beamfd)
26	
27	# --- Tests ------------------------------------------------------------------
28	enable_testing()
29	add_subdirectory(tests)
30

ls: cannot access 'build': No such file or directory
---

[stdout]
ls: cannot access 'build': No such file or directory
---

/opt/beamfd/build/libbeamfd.a

[stdout]
/opt/beamfd/build/libbeamfd.a

pinned-pinned UDL                  err(N)=1.450e-06 err(2N)=3.625e-07 ratio=4.00 order=2.00
cantilever tip load                err(N)=2.813e-06 err(2N)=7.031e-07 ratio=4.00 order=2.00
cantilever tip moment              err(N)=2.006e-13 err(2N)=3.709e-12 ratio=0.05 order=-4.21
clamped-clamped UDL                err(N)=8.121e-07 err(2N)=2.070e-07 ratio=3.92 order=1.97
variable-EI + foundation (MMS)     err(N)=1.767e-04 err(2N)=5.323e-05 ratio=3.32 order=1.73
free-free on foundation (MMS)      err(N)=4.473e+00 err(2N)=2.131e+00 ratio=2.10 order=1.07

[stdout]
pinned-pinned UDL                  err(N)=1.450e-06 err(2N)=3.625e-07 ratio=4.00 order=2.00
cantilever tip load                err(N)=2.813e-06 err(2N)=7.031e-07 ratio=4.00 order=2.00
cantilever tip moment              err(N)=2.006e-13 err(2N)=3.709e-12 ratio=0.05 order=-4.21
clamped-clamped UDL                err(N)=8.121e-07 err(2N)=2.070e-07 ratio=3.92 order=1.97
variable-EI + foundation (MMS)     err(N)=1.767e-04 err(2N)=5.323e-05 ratio=3.32 order=1.73
free-free on foundation (MMS)      err(N)=4.473e+00 err(2N)=2.131e+00 ratio=2.10 order=1.07

BUILD_OK
pinned-pinned UDL                  err(N)=1.450e-06 err(2N)=3.625e-07 ratio=4.00 order=2.00
cantilever tip load                err(N)=2.813e-06 err(2N)=7.031e-07 ratio=4.00 order=2.00
cantilever tip moment              err(N)=1.987e-13 err(2N)=3.700e-12 ratio=0.05 order=-4.22
clamped-clamped UDL                err(N)=8.121e-07 err(2N)=2.070e-07 ratio=3.92 order=1.97
variable-EI + foundation (MMS)     err(N)=1.767e-04 err(2N)=5.323e-05 ratio=3.32 order=1.73
variable-EI + foundation (MMS) fine err(N)=5.323e-05 err(2N)=1.444e-05 ratio=3.69 order=1.88
free-free on foundation (MMS)      err(N)=6.320e-02 err(2N)=1.709e-02 ratio=3.70 order=1.89
free-free on foundation (MMS) fine err(N)=1.709e-02 err(2N)=4.435e-03 ratio=3.85 order=1.95
clamped-free curved (MMS)          err(N)=9.103e-04 err(2N)=2.522e-04 ratio=3.61 order=1.85
clamped-free curved (MMS) fine     err(N)=2.522e-04 err(2N)=6.612e-05 ratio=3.81 order=1.93

[stdout]
BUILD_OK
pinned-pinned UDL                  err(N)=1.450e-06 err(2N)=3.625e-07 ratio=4.00 order=2.00
cantilever tip load                err(N)=2.813e-06 err(2N)=7.031e-07 ratio=4.00 order=2.00
cantilever tip moment              err(N)=1.987e-13 err(2N)=3.700e-12 ratio=0.05 order=-4.22
clamped-clamped UDL                err(N)=8.121e-07 err(2N)=2.070e-07 ratio=3.92 order=1.97
variable-EI + foundation (MMS)     err(N)=1.767e-04 err(2N)=5.323e-05 ratio=3.32 order=1.73
variable-EI + foundation (MMS) fine err(N)=5.323e-05 err(2N)=1.444e-05 ratio=3.69 order=1.88
free-free on foundation (MMS)      err(N)=6.320e-02 err(2N)=1.709e-02 ratio=3.70 order=1.89
free-free on foundation (MMS) fine err(N)=1.709e-02 err(2N)=4.435e-03 ratio=3.85 order=1.95
clamped-free curved (MMS)          err(N)=9.103e-04 err(2N)=2.522e-04 ratio=3.61 order=1.85
clamped-free curved (MMS) fine     err(N)=2.522e-04 err(2N)=6.612e-05 ratio=3.81 order=1.93

NO WARNINGS
---CTEST---
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]
NO WARNINGS
---CTEST---
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 2
cleaned up scratch files
ls: cannot access 'verify.cpp': No such file or directory

[error] tool reported failure

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	    // Uniform grid x_i = i*h, i = 0..N, with n = N+1 nodal unknowns w_i.
41	    //
42	    // Governing equation (conservative / divergence form):
43	    //
44	    //     d^2/dx^2 ( EI(x) w''(x) ) + k(x) w(x) = q(x)
45	    //
46	    // Discretization (second-order accurate):
47	    //   * Interior rows i = 2..N-2 enforce the conservative bending operator as
48	    //     the second difference of the nodal bending moment m_j = EI_j w''_j:
49	    //
50	    //         L_i = ( m_{i-1} - 2 m_i + m_{i+1} ) / h^2 ,
51	    //         m_j = EI_j ( w_{j-1} - 2 w_j + w_{j+1} ) / h^2 ,
52	    //
53	    //     giving a five-point stencil that reduces to the usual biharmonic
54	    //     stencil EI/h^4 * [1 -4 6 -4 1] when EI is constant. The discrete
55	    //     equation is L_i + k_i w_i = q_i.
56	    //   * The first/last two rows are the boundary conditions, written as
57	    //     one-sided difference stencils consistent to the same order. The
58	    //     natural (free-end) rows reuse exactly the discrete M and V operators
59	    //     of the recovery helpers, so their RHS are end_moment / applied_end_shear.
60	    const std::size_t n = beam_.num_nodes();
61	    const std::size_t N = beam_.num_intervals();  // n - 1
62	    const double h = beam_.dx();
63	    const double h2 = h * h;
64	    const double h4 = h2 * h2;
65	
66	    BeamSystem sys;
67	    sys.A = DenseMatrix(n, n);
68	    sys.b.assign(n, 0.0);
69	    DenseMatrix& A = sys.A;
70	    std::vector<double>& b = sys.b;
71	
72	    // ---- Interior conservative bending rows: i = 2 .. N-2 ------------------
73	    for (std::size_t i = 2; i + 2 <= N; ++i) {
74	        const double EIm = beam_.EI_at(i - 1);
75	        const double EIc = beam_.EI_at(i);
76	        const double EIp = beam_.EI_at(i + 1);
77	        A(i, i - 2) += EIm / h4;
78	        A(i, i - 1) += (-2.0 * EIm - 2.0 * EIc) / h4;
79	        A(i, i) += (EIm + 4.0 * EIc + EIp) / h4 + beam_.k_at(i);
80	        A(i, i + 1) += (-2.0 * EIc - 2.0 * EIp) / h4;
81	        A(i, i + 2) += EIp / h4;
82	        b[i] = beam_.q_at(i);
83	    }
84	
85	    // ---- Left boundary conditions: rows 0 and 1 ---------------------------
86	    switch (beam_.left()) {
87	        case Support::Clamped: {
88	            // w(0) = 0
89	            A(0, 0) = 1.0;
90	            b[0] = 0.0;
91	            // w'(0) = 0 via second-order forward difference.
92	            A(1, 0) = -3.0 / (2.0 * h);
93	            A(1, 1) = 4.0 / (2.0 * h);
94	            A(1, 2) = -1.0 / (2.0 * h);
95	            b[1] = 0.0;
96	            break;
97	        }
98	        case Support::Pinned: {
99	            // w(0) = 0
100	            A(0, 0) = 1.0;
101	            b[0] = 0.0;
102	            // M(0) = EI w''(0) = 0; one-sided second derivative (scaled to O(1)).
103	            A(1, 0) = 2.0;
104	            A(1, 1) = -5.0;
105	            A(1, 2) = 4.0;
106	            A(1, 3) = -1.0;
107	            b[1] = 0.0;
108	            break;
109	        }
110	        case Support::Free: {
111	            const double EI0 = beam_.EI_at(0);
112	            const double EI1 = beam_.EI_at(1);
113	            const double EI2 = beam_.EI_at(2);
114	            // Row 0: M(0) = EI0 * (2 w0 - 5 w1 + 4 w2 - w3)/h^2 = end_moment.
115	            A(0, 0) = EI0 * 2.0 / h2;
116	            A(0, 1) = EI0 * -5.0 / h2;
117	            A(0, 2) = EI0 * 4.0 / h2;
118	            A(0, 3) = EI0 * -1.0 / h2;
119	            b[0] = beam_.end_moment(true);
120	            // Row 1: V(0) = (EI w'')'(0) = EI'(0) w''(0) + EI(0) w'''(0), each
121	            // factor a second-order one-sided difference (product rule keeps the
122	            // whole shear operator O(h^2), unlike a composite nodal-moment form):
123	            //   EI'(0) ~ (-3 EI0 + 4 EI1 - EI2)/(2h)
124	            //   w''(0) ~ (2 w0 - 5 w1 + 4 w2 - w3)/h^2
125	            //   w'''(0) ~ (-5 w0 + 18 w1 - 24 w2 + 14 w3 - 3 w4)/(2 h^3)
126	            const double dEI0 = (-3.0 * EI0 + 4.0 * EI1 - EI2) / (2.0 * h);
127	            const double cm = dEI0 / h2;                 // EI'(0) * w''(0) stencil
128	            const double cs = EI0 / (2.0 * h * h2);      // EI(0) * w'''(0) stencil
129	            A(1, 0) += cm * 2.0 + cs * -5.0;
130	            A(1, 1) += cm * -5.0 + cs * 18.0;
131	            A(1, 2) += cm * 4.0 + cs * -24.0;
132	            A(1, 3) += cm * -1.0 + cs * 14.0;
133	            A(1, 4) += cs * -3.0;
134	            b[1] = beam_.applied_end_shear(true);
135	            break;
136	        }
137	    }
138	
139	    // ---- Right boundary conditions: rows N-1 and N ------------------------
140	    switch (beam_.right()) {
141	        case Support::Clamped: {
142	            // w(L) = 0
143	            A(N, N) = 1.0;
144	            b[N] = 0.0;
145	            // w'(L) = 0 via second-order backward difference.
146	            A(N - 1, N) = 3.0 / (2.0 * h);
147	            A(N - 1, N - 1) = -4.0 / (2.0 * h);
148	            A(N - 1, N - 2) = 1.0 / (2.0 * h);
149	            b[N - 1] = 0.0;
150	            break;
151	        }
152	        case Support::Pinned: {
153	            // w(L) = 0
154	            A(N, N) = 1.0;
155	            b[N] = 0.0;
156	            // M(L) = EI w''(L) = 0; one-sided second derivative (scaled to O(1)).
157	            A(N - 1, N) = 2.0;
158	            A(N - 1, N - 1) = -5.0;
159	            A(N - 1, N - 2) = 4.0;
160	            A(N - 1, N - 3) = -1.0;
161	            b[N - 1] = 0.0;
162	            break;
163	        }
164	        case Support::Free: {
165	            const double EIN = beam_.EI_at(N);
166	            const double EINm1 = beam_.EI_at(N - 1);
167	            const double EINm2 = beam_.EI_at(N - 2);
168	            // Row N: M(L) = EIN * (2 wN - 5 w_{N-1} + 4 w_{N-2} - w_{N-3})/h^2.
169	            A(N, N) = EIN * 2.0 / h2;
170	            A(N, N - 1) = EIN * -5.0 / h2;
171	            A(N, N - 2) = EIN * 4.0 / h2;
172	            A(N, N - 3) = EIN * -1.0 / h2;
173	            b[N] = beam_.end_moment(false);
174	            // Row N-1: V(L) = EI'(L) w''(L) + EI(L) w'''(L), second-order
175	            // one-sided (backward) differences (mirror of the left free end):
176	            //   EI'(L) ~ (3 EIN - 4 EI_{N-1} + EI_{N-2})/(2h)
177	            //   w''(L) ~ (2 wN - 5 w_{N-1} + 4 w_{N-2} - w_{N-3})/h^2
178	            //   w'''(L) ~ (5 wN - 18 w_{N-1} + 24 w_{N-2} - 14 w_{N-3} + 3 w_{N-4})/(2 h^3)
179	            const double dEIN = (3.0 * EIN - 4.0 * EINm1 + EINm2) / (2.0 * h);
180	            const double cm = dEIN / h2;             // EI'(L) * w''(L) stencil
181	            const double cs = EIN / (2.0 * h * h2);  // EI(L) * w'''(L) stencil
182	            A(N - 1, N) += cm * 2.0 + cs * 5.0;
183	            A(N - 1, N - 1) += cm * -5.0 + cs * -18.0;
184	            A(N - 1, N - 2) += cm * 4.0 + cs * 24.0;
185	            A(N - 1, N - 3) += cm * -1.0 + cs * -14.0;
186	            A(N - 1, N - 4) += cs * 3.0;
187	            b[N - 1] = beam_.applied_end_shear(false);
188	            break;
189	        }
190	    }
191	
192	    return sys;
193	}
194	
195	BeamResult BeamSolver::solve() const {
196	    // Validate the support/foundation configuration first; let mechanism errors
197	    // propagate. Then assemble and solve A w = b (a singular system surfaces as
198	    // an exception from solve_general, which also propagates).
199	    beam_.validate();
200	    const BeamSystem sys = assemble_system();

BUILD CLEAN
100% tests passed, 0 tests failed out of 1

Total Test time (real) =   0.00 sec

[stdout]
BUILD CLEAN
100% tests passed, 0 tests failed out of 1

Total Test time (real) =   0.00 sec