tests/test.sh exit code (0 → resolved); the classification below is post-hoc and cannot change it./opt/beamfd/include/beamfd/solver.hpp
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
/opt/beamfd/src/solver.cpp
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
/opt/beamfd/include/beamfd/beam.hpp
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
/opt/beamfd/include/beamfd/linalg.hpp
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
/opt/beamfd/include/beamfd/types.hpp
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
/opt/beamfd/src/beam.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
**/*.cpp
tests/test_smoke.cpp
**/*
CMakeLists.txt examples/cantilever_tip.beam include/beamfd/linalg.hpp src/linalg.cpp src/main.cpp tests/test_util.hpp include/beamfd/types.hpp include/beamfd/beam.hpp src/beam.cpp include/beamfd/solver.hpp src/solver.cpp src/io.cpp include/beamfd/io.hpp tests/test_smoke.cpp tests/CMakeLists.txt
/opt/beamfd/tests/test_smoke.cpp
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
/opt/beamfd/tests/test_util.hpp
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
/opt/beamfd/examples/cantilever_tip.beam
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
/opt/beamfd/src/linalg.cpp
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
/opt/beamfd/src/solver.cpp
BeamSystem BeamSolver::assemble_system() const {
// TODO(candidate): build and return the finite-difference system A w = b.
const std::size_t n = beam_.num_nodes();
BeamSystem sys;
sys.A = DenseMatrix(n, n); // all zeros -> singular / wrong
sys.b.assign(n, 0.0);
return sys;
}
BeamResult BeamSolver::solve() const {
// TODO(candidate): validate, assemble, solve, return the deflection field.
const std::size_t n = beam_.num_nodes();
BeamResult result;
result.deflection.assign(n, 0.0); // trivially zero -> wrong
return result;
}// ----------------------------------------------------------------------------
// Finite-difference discretization of
//
// d^2/dx^2 ( EI(x) w''(x) ) + k(x) w(x) = q(x)
//
// on the uniform grid x_i = i*dx, i = 0..N (n = N+1 nodes), second-order
// accurate throughout.
//
// The fourth-order operator is built by nesting two central second differences:
// the nodal moment is M_j = EI_j * w''_j with the *same* second-difference
// stencils used by the recovery helpers (central in the interior, a 2nd-order
// one-sided form at an end), and the governing equation is the second
// difference of that moment field,
//
// (M_{i-1} - 2 M_i + M_{i+1}) / dx^2 + k_i w_i = q_i,
//
// expanded into a five-point [w_{i-2} .. w_{i+2}] stencil with variable EI:
//
// w_{i-2}: EI_{i-1} \
// w_{i-1}: -2(EI_{i-1} + EI_i) | all /dx^4
// w_i : EI_{i-1} + 4 EI_i + EI_{i+1} |
// w_{i+1}: -2(EI_i + EI_{i+1}) |
// w_{i+2}: EI_{i+1} /
//
// This interior equation is enforced at nodes i = 2 .. N-2. The remaining four
// rows (0, 1, N-1, N) carry the two boundary conditions at each end:
//
// * "outer" row (0 / N): the deflection-or-moment condition
// clamped / pinned : w = 0
// free : EI*w'' = M_applied
// * "inner" row (1 / N-1): the slope-or-shear condition
// clamped : w' = 0
// pinned : EI*w'' = 0 (zero moment)
// free : (EI*w'')' = V_applied (shear)
//
// The slope, moment and shear stencils are all 2nd-order one-sided forms
// consistent with the interior scheme and with the recovery helpers.
BeamSystem BeamSolver::assemble_system() const {
const std::size_t n = beam_.num_nodes();
const std::size_t N = beam_.num_intervals(); // last node index
const double h = beam_.dx();
const double h2 = h * h;
const double h4 = h2 * h2;
BeamSystem sys;
sys.A = DenseMatrix(n, n);
sys.b.assign(n, 0.0);
DenseMatrix& A = sys.A;
std::vector<double>& b = sys.b;
// ---- Interior governing equation at nodes 2 .. N-2 ---------------------
for (std::size_t i = 2; i + 2 <= N; ++i) {
const double eim = beam_.EI_at(i - 1);
const double eii = beam_.EI_at(i);
const double eip = beam_.EI_at(i + 1);
A(i, i - 2) += eim / h4;
A(i, i - 1) += -2.0 * (eim + eii) / h4;
A(i, i) += (eim + 4.0 * eii + eip) / h4 + beam_.k_at(i);
A(i, i + 1) += -2.0 * (eii + eip) / h4;
A(i, i + 2) += eip / h4;
b[i] = beam_.q_at(i);
}
// ---- Boundary rows -----------------------------------------------------
// The four stencils below are written for the left end in terms of the
// four boundary nodes (e0, e1, e2, e3) = (0, 1, 2, 3); the right end uses
// the mirror image (N, N-1, N-2, N-3) with the sign of the (odd-order)
// shear stencil flipped.
auto fill_end = [&](bool left_end) {
const Support sup = left_end ? beam_.left() : beam_.right();
// Node indices walking inward from the end, and the row indices for
// the outer (deflection/moment) and inner (slope/shear) conditions.
const std::size_t e0 = left_end ? 0 : N;
const std::size_t e1 = left_end ? 1 : N - 1;
const std::size_t e2 = left_end ? 2 : N - 2;
const std::size_t e3 = left_end ? 3 : N - 3;
const std::size_t outer_row = e0;
const std::size_t inner_row = e1;
const double ei0 = beam_.EI_at(e0);
const double ei1 = beam_.EI_at(e1);
const double ei2 = beam_.EI_at(e2);
// Forward (+1) at the left end, backward (-1) at the right end, used to
// orient the one-sided first-derivative (slope/shear) stencils.
const double dir = left_end ? 1.0 : -1.0;
// Second-derivative (w'') one-sided stencil at the end node:
// (2 w_e0 - 5 w_e1 + 4 w_e2 - w_e3) / h^2 (2nd order)
auto set_w2 = [&](std::size_t row, double scale, double rhs) {
A(row, e0) += scale * 2.0 / h2;
A(row, e1) += scale * -5.0 / h2;
A(row, e2) += scale * 4.0 / h2;
A(row, e3) += scale * -1.0 / h2;
b[row] = rhs;
};
// --- Outer row: deflection (essential) or moment (natural) ---------
if (sup == Support::Free) {
// EI * w'' = M_applied
set_w2(outer_row, ei0, beam_.end_moment(left_end));
} else {
// w = 0
A(outer_row, e0) += 1.0;
b[outer_row] = 0.0;
}
// --- Inner row: slope (clamped) / moment (pinned) / shear (free) ---
if (sup == Support::Clamped) {
// w' = 0, 2nd-order one-sided: (-3 w_e0 + 4 w_e1 - w_e2)/(2h) = 0
A(inner_row, e0) += dir * -3.0 / (2.0 * h);
A(inner_row, e1) += dir * 4.0 / (2.0 * h);
A(inner_row, e2) += dir * -1.0 / (2.0 * h);
b[inner_row] = 0.0;
} else if (sup == Support::Pinned) {
// EI * w'' = 0 (zero moment)
set_w2(inner_row, ei0, 0.0);
} else { // Free: (EI w'')' = V_applied
// Shear = M'(end) via a one-sided first difference of the nodal
// moment field M_j = EI_j w''_j:
// M_e0 = EI_e0 (2 w_e0 - 5 w_e1 + 4 w_e2 - w_e3)/h^2 (one-sided)
// M_e1 = EI_e1 ( w_e0 - 2 w_e1 + w_e2 )/h^2 (central)
// M_e2 = EI_e2 ( w_e1 - 2 w_e2 + w_e3 )/h^2 (central)
// M'(end) ~ (-3 M_e0 + 4 M_e1 - M_e2)/(2h)
// Combined coefficients (factor 1/(2 h^3)):
const double f = dir / (2.0 * h * h2);
A(inner_row, e0) += f * (-6.0 * ei0 + 4.0 * ei1);
A(inner_row, e1) += f * (15.0 * ei0 - 8.0 * ei1 - ei2);
A(inner_row, e2) += f * (-12.0 * ei0 + 4.0 * ei1 + 2.0 * ei2);
A(inner_row, e3) += f * (3.0 * ei0 - ei2);
b[inner_row] = beam_.applied_end_shear(left_end);
}
};
fill_end(true);
fill_end(false);
return sys;
}
BeamResult BeamSolver::solve() const {
beam_.validate(); // propagate mechanism / ill-posed configuration errors
BeamSystem sys = assemble_system();
BeamResult result;
result.deflection = solve_general(sys.A, sys.b); // propagates if singular
return result;
}Configure, build, and run tests
cmake -S . -B build -DCMAKE_BUILD_TYPE=Release > /tmp/cmake.log 2>&1 && cmake --build build --parallel > /tmp/build.log 2>&1 && echo BUILD_OK && cd build && ctest --output-on-failure 2>&1 | tail -20
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 secVerify convergence against closed-form solutions
cat > /tmp/verify.cpp <<'EOF'
#include <cmath>
#include <cstdio>
#include <vector>
#include "beamfd/beam.hpp"
#include "beamfd/solver.hpp"
using namespace beamfd;
double maxabs_err(const std::vector<double>&w, double L, std::size_t N, double(*ex)(double,double), double L_) {
double e=0; for(std::size_t i=0;i<=N;++i){double x=L*i/N; e=std::max(e,std::fabs(w[i]-ex(x,L_)));} return e;
}
// exact solutions
double ss_udl(double x,double L){ /*EI=1,q=1*/ return (x*(L*L*L - 2*L*x*x + x*x*x))/24.0; }
double cc_udl(double x,double L){ return (x*x*(L-x)*(L-x))/24.0; }
double cant_tip(double x,double L){ /*P=1 down, EI=1: w=P/(6EI)(3L x^2 - x^3)*/ return (3*L*x*x - x*x*x)/6.0; }
void conv(const char*name, double(*build_err)(std::size_t)){
printf("%-22s ", name);
double prev=0; for(int k=0;k<4;++k){std::size_t N=20*(1<<k); double e=build_err(N);
printf("N=%4zu err=%.3e", N, e); if(k){printf(" rate=%.2f",std::log(prev/e)/std::log(2.0));} printf(" "); prev=e;}
printf("\n");
}
double err_ss(std::size_t N){ Beam b(3.0,1.0,N,Support::Pinned,Support::Pinned); b.set_distributed_load(1.0);
auto w=BeamSolver(b).solve().deflection; double e=0;for(std::size_t i=0;i<=N;++i){double x=3.0*i/N;e=std::max(e,std::fabs(w[i]-ss_udl(x,3.0)));}return e;}
double err_cc(std::size_t N){ Beam b(3.0,1.0,N,Support::Clamped,Support::Clamped); b.set_distributed_load(1.0);
auto w=BeamSolver(b).solve().deflection; double e=0;for(std::size_t i=0;i<=N;++i){double x=3.0*i/N;e=std::max(e,std::fabs(w[i]-cc_udl(x,3.0)));}return e;}
double err_cant(std::size_t N){ Beam b(3.0,1.0,N,Support::Clamped,Support::Free); b.set_end_load(1.0,false);
auto w=BeamSolver(b).solve().deflection; double e=0;for(std::size_t i=0;i<=N;++i){double x=3.0*i/N;e=std::max(e,std::fabs(w[i]-cant_tip(x,3.0)));}return e;}
int main(){
conv("pinned-pinned UDL", err_ss);
conv("clamped-clamped UDL", err_cc);
conv("cantilever tip load", err_cant);
return 0;
}
EOF
g++ -std=c++17 -I/opt/beamfd/include /tmp/verify.cpp /opt/beamfd/src/beam.cpp /opt/beamfd/src/solver.cpp /opt/beamfd/src/linalg.cpp -o /tmp/verify && /tmp/verifypinned-pinned UDL N= 20 err=2.320e-02 N= 40 err=5.801e-03 rate=2.00 N= 80 err=1.450e-03 rate=2.00 N= 160 err=3.625e-04 rate=2.00 clamped-clamped UDL N= 20 err=7.805e-03 N= 40 err=2.030e-03 rate=1.94 N= 80 err=5.175e-04 rate=1.97 N= 160 err=1.306e-04 rate=1.99 cantilever tip load N= 20 err=2.250e-02 N= 40 err=5.625e-03 rate=2.00 N= 80 err=1.406e-03 rate=2.00 N= 160 err=3.515e-04 rate=2.00 [stdout] pinned-pinned UDL N= 20 err=2.320e-02 N= 40 err=5.801e-03 rate=2.00 N= 80 err=1.450e-03 rate=2.00 N= 160 err=3.625e-04 rate=2.00 clamped-clamped UDL N= 20 err=7.805e-03 N= 40 err=2.030e-03 rate=1.94 N= 80 err=5.175e-04 rate=1.97 N= 160 err=1.306e-04 rate=1.99 cantilever tip load N= 20 err=2.250e-02 N= 40 err=5.625e-03 rate=2.00 N= 80 err=1.406e-03 rate=2.00 N= 160 err=3.515e-04 rate=2.00
Verify variable-EI and foundation manufactured solutions
cat > /tmp/verify2.cpp <<'EOF'
#include <cmath>
#include <cstdio>
#include <vector>
#include "beamfd/beam.hpp"
#include "beamfd/solver.hpp"
using namespace beamfd;
// Manufactured solution w = x^2 (L-x)^2 ; EI = p+qx ; k = c+d x.
// w'' = 2L^2 -12Lx +12x^2 ; M=EI w'' ; M'' = 24p - 24qL + 72 q x ; load = M'' + k w.
struct MS { double L,p,q,c,d;
double w(double x)const{return x*x*(L-x)*(L-x);}
double w2(double x)const{return 2*L*L-12*L*x+12*x*x;}
double EI(double x)const{return p+q*x;}
double k(double x)const{return c+d*x;}
double M(double x)const{return EI(x)*w2(x);}
double Mpp(double x)const{return 24*p-24*q*L+72*q*x;}
double Mp(double x)const{double g=w2(x),gp=-12*L+24*x;return q*g+EI(x)*gp;}
double load(double x)const{return Mpp(x)+k(x)*w(x);}
};
double run(const MS&ms, std::size_t N, bool clampedFree){
double L=ms.L;
Beam b(L,1.0,N, Support::Clamped, clampedFree?Support::Free:Support::Clamped);
std::vector<double> ei(N+1),kk(N+1),q(N+1);
for(std::size_t i=0;i<=N;++i){double x=L*i/N; ei[i]=ms.EI(x); kk[i]=ms.k(x); q[i]=ms.load(x);}
b.set_ei_nodal(ei); b.set_foundation_nodal(kk); b.set_q_nodal(q);
if(clampedFree){
// free at right: EI w'' = end_moment ; (EI w'')' = applied_end_shear = -force
b.set_end_moment(ms.M(L), false);
b.set_end_load(-ms.Mp(L), false);
}
auto w=BeamSolver(b).solve().deflection;
double e=0;for(std::size_t i=0;i<=N;++i){double x=L*i/N;e=std::max(e,std::fabs(w[i]-ms.w(x)));}
return e;
}
void conv(const char*name, const MS&ms, bool cf){
printf("%-28s ",name); double prev=0;
for(int k=0;k<4;++k){std::size_t N=20*(1<<k); double e=run(ms,N,cf);
printf("err=%.3e",e); if(k)printf("(r=%.2f)",std::log(prev/e)/std::log(2.0)); printf(" "); prev=e;}
printf("\n");
}
int main(){
MS varEI{2.0, 1000.0, 500.0, 0.0, 0.0}; // linear EI, no foundation
MS withFound{2.0, 1000.0, 500.0, 300.0, 200.0}; // varEI + linear foundation
MS foundOnly{2.0, 800.0, 0.0, 400.0, 0.0};
conv("varEI clamped-clamped", varEI, false);
conv("varEI+found clamp-clamp", withFound, false);
conv("varEI clamped-FREE", varEI, true);
conv("varEI+found clamped-FREE", withFound, true);
conv("const EI+found c-c", foundOnly, false);
return 0;
}
EOF
g++ -std=c++17 -I/opt/beamfd/include /tmp/verify2.cpp /opt/beamfd/src/beam.cpp /opt/beamfd/src/solver.cpp /opt/beamfd/src/linalg.cpp -o /tmp/verify2 && /tmp/verify2varEI clamped-clamped err=3.700e-02 err=9.625e-03(r=1.94) err=2.453e-03(r=1.97) err=6.191e-04(r=1.99) varEI+found clamp-clamp err=3.657e-02 err=9.510e-03(r=1.94) err=2.424e-03(r=1.97) err=6.117e-04(r=1.99) varEI clamped-FREE err=1.419e+01 err=7.479e+00(r=0.92) err=3.833e+00(r=0.96) err=1.940e+00(r=0.98) varEI+found clamped-FREE err=9.730e+00 err=4.811e+00(r=1.02) err=2.386e+00(r=1.01) err=1.187e+00(r=1.01) const EI+found c-c err=3.636e-02 err=9.455e-03(r=1.94) err=2.410e-03(r=1.97) err=6.081e-04(r=1.99) [stdout] varEI clamped-clamped err=3.700e-02 err=9.625e-03(r=1.94) err=2.453e-03(r=1.97) err=6.191e-04(r=1.99) varEI+found clamp-clamp err=3.657e-02 err=9.510e-03(r=1.94) err=2.424e-03(r=1.97) err=6.117e-04(r=1.99) varEI clamped-FREE err=1.419e+01 err=7.479e+00(r=0.92) err=3.833e+00(r=0.96) err=1.940e+00(r=0.98) varEI+found clamped-FREE err=9.730e+00 err=4.811e+00(r=1.02) err=2.386e+00(r=1.01) err=1.187e+00(r=1.01) const EI+found c-c err=3.636e-02 err=9.455e-03(r=1.94) err=2.410e-03(r=1.97) err=6.081e-04(r=1.99)
Test constant-EI clamped-free manufactured solution
cat >> /tmp/verify2.cpp <<'EOF'
EOF
# add a const-EI clamped-free case by editing main via a quick second program
cat > /tmp/verify3.cpp <<'EOF'
#include <cmath>
#include <cstdio>
#include <vector>
#include "beamfd/beam.hpp"
#include "beamfd/solver.hpp"
using namespace beamfd;
struct MS { double L,p,q,c,d;
double w(double x)const{return x*x*(L-x)*(L-x);}
double w2(double x)const{return 2*L*L-12*L*x+12*x*x;}
double EI(double x)const{return p+q*x;}
double k(double x)const{return c+d*x;}
double M(double x)const{return EI(x)*w2(x);}
double Mpp(double x)const{return 24*p-24*q*L+72*q*x;}
double Mp(double x)const{double g=w2(x),gp=-12*L+24*x;return q*g+EI(x)*gp;}
double load(double x)const{return Mpp(x)+k(x)*w(x);}
};
double run(const MS&ms, std::size_t N){
double L=ms.L; Beam b(L,1.0,N, Support::Clamped, Support::Free);
std::vector<double> ei(N+1),kk(N+1),q(N+1);
for(std::size_t i=0;i<=N;++i){double x=L*i/N; ei[i]=ms.EI(x); kk[i]=ms.k(x); q[i]=ms.load(x);}
b.set_ei_nodal(ei); b.set_foundation_nodal(kk); b.set_q_nodal(q);
b.set_end_moment(ms.M(L), false); b.set_end_load(-ms.Mp(L), false);
auto w=BeamSolver(b).solve().deflection; double e=0,wm=0;
for(std::size_t i=0;i<=N;++i){double x=L*i/N;e=std::max(e,std::fabs(w[i]-ms.w(x)));wm=std::max(wm,std::fabs(ms.w(x)));}
return e;
}
void conv(const char*name,const MS&ms){printf("%-26s ",name);double prev=0;
for(int k=0;k<4;++k){std::size_t N=20*(1<<k);double e=run(ms,N);printf("err=%.3e",e);
if(k)printf("(r=%.2f)",std::log(prev/e)/std::log(2.0));printf(" ");prev=e;}printf("\n");}
int main(){
conv("constEI clamped-FREE", MS{2.0,1000.0,0.0,0.0,0.0});
return 0;
}
EOF
g++ -std=c++17 -I/opt/beamfd/include /tmp/verify3.cpp /opt/beamfd/src/beam.cpp /opt/beamfd/src/solver.cpp /opt/beamfd/src/linalg.cpp -o /tmp/verify3 && /tmp/verify3constEI clamped-FREE err=9.284e+00 err=4.726e+00(r=0.97) err=2.382e+00(r=0.99) err=1.196e+00(r=0.99) [stdout] constEI clamped-FREE err=9.284e+00 err=4.726e+00(r=0.97) err=2.382e+00(r=0.99) err=1.196e+00(r=0.99)
/opt/beamfd/src/solver.cpp
39 // ----------------------------------------------------------------------------
40 // Finite-difference discretization of
41 //
42 // d^2/dx^2 ( EI(x) w''(x) ) + k(x) w(x) = q(x)
43 //
44 // on the uniform grid x_i = i*dx, i = 0..N (n = N+1 nodes), second-order
45 // accurate throughout.
46 //
47 // The fourth-order operator is built by nesting two central second differences:
48 // the nodal moment is M_j = EI_j * w''_j with the *same* second-difference
49 // stencils used by the recovery helpers (central in the interior, a 2nd-order
50 // one-sided form at an end), and the governing equation is the second
51 // difference of that moment field,
52 //
53 // (M_{i-1} - 2 M_i + M_{i+1}) / dx^2 + k_i w_i = q_i,
54 //
55 // expanded into a five-point [w_{i-2} .. w_{i+2}] stencil with variable EI:
56 //
57 // w_{i-2}: EI_{i-1} \
58 // w_{i-1}: -2(EI_{i-1} + EI_i) | all /dx^4
59 // w_i : EI_{i-1} + 4 EI_i + EI_{i+1} |
60 // w_{i+1}: -2(EI_i + EI_{i+1}) |
61 // w_{i+2}: EI_{i+1} /
62 //
63 // This interior equation is enforced at nodes i = 2 .. N-2. The remaining four
64 // rows (0, 1, N-1, N) carry the two boundary conditions at each end:
65 //
66 // * "outer" row (0 / N): the deflection-or-moment condition
67 // clamped / pinned : w = 0
68 // free : EI*w'' = M_applied
69 // * "inner" row (1 / N-1): the slope-or-shear condition
70 // clamped : w' = 0
71 // pinned : EI*w'' = 0 (zero moment)
72 // free : (EI*w'')' = V_applied (shear)
73 //
74 // The slope, moment and shear stencils are all 2nd-order one-sided forms
75 // consistent with the interior scheme and with the recovery helpers.
76
77 BeamSystem BeamSolver::assemble_system() const {
78 const std::size_t n = beam_.num_nodes();
79 const std::size_t N = beam_.num_intervals(); // last node index
80 const double h = beam_.dx();
81 const double h2 = h * h;
82 const double h4 = h2 * h2;
83
84 BeamSystem sys;
85 sys.A = DenseMatrix(n, n);
86 sys.b.assign(n, 0.0);
87 DenseMatrix& A = sys.A;
88 std::vector<double>& b = sys.b;
89
90 // ---- Interior governing equation at nodes 2 .. N-2 ---------------------
91 for (std::size_t i = 2; i + 2 <= N; ++i) {
92 const double eim = beam_.EI_at(i - 1);
93 const double eii = beam_.EI_at(i);
/opt/beamfd/src/solver.cpp
39 // ----------------------------------------------------------------------------
40 // Finite-difference discretization of
41 //
42 // d^2/dx^2 ( EI(x) w''(x) ) + k(x) w(x) = q(x)
43 //
44 // on the uniform grid x_i = i*dx, i = 0..N (n = N+1 nodes), second-order
45 // accurate throughout.
46 //
47 // The fourth-order operator is built by nesting two central second differences:
48 // the nodal moment is M_j = EI_j * w''_j with the *same* second-difference
49 // stencils used by the recovery helpers (central in the interior, a 2nd-order
50 // one-sided form at an end), and the governing equation is the second
51 // difference of that moment field,
52 //
53 // (M_{i-1} - 2 M_i + M_{i+1}) / dx^2 + k_i w_i = q_i,
54 //
55 // expanded into a five-point [w_{i-2} .. w_{i+2}] stencil with variable EI:
56 //
57 // w_{i-2}: EI_{i-1} \
58 // w_{i-1}: -2(EI_{i-1} + EI_i) | all /dx^4
59 // w_i : EI_{i-1} + 4 EI_i + EI_{i+1} |
60 // w_{i+1}: -2(EI_i + EI_{i+1}) |
61 // w_{i+2}: EI_{i+1} /
62 //
63 // This interior equation is enforced at nodes i = 2 .. N-2. The remaining four
64 // rows (0, 1, N-1, N) carry the two boundary conditions at each end:
65 //
66 // * "outer" row (0 / N): the deflection-or-moment condition
67 // clamped / pinned : w = 0
68 // free : EI*w'' = M_applied
69 // * "inner" row (1 / N-1): the slope-or-shear condition
70 // clamped : w' = 0
71 // pinned : EI*w'' = 0 (zero moment)
72 // free : (EI*w'')' = V_applied (shear)
73 //
74 // The slope, moment and shear stencils are all 2nd-order one-sided forms
75 // consistent with the interior scheme and with the recovery helpers.
76
77 BeamSystem BeamSolver::assemble_system() const {
78 const std::size_t n = beam_.num_nodes();
79 const std::size_t N = beam_.num_intervals(); // last node index
80 const double h = beam_.dx();
81 const double h2 = h * h;
82 const double h4 = h2 * h2;
83
84 BeamSystem sys;
85 sys.A = DenseMatrix(n, n);
86 sys.b.assign(n, 0.0);
87 DenseMatrix& A = sys.A;
88 std::vector<double>& b = sys.b;
89
90 // ---- Interior governing equation at nodes 2 .. N-2 ---------------------
91 for (std::size_t i = 2; i + 2 <= N; ++i) {
92 const double eim = beam_.EI_at(i - 1);
93 const double eii = beam_.EI_at(i);
94 const double eip = beam_.EI_at(i + 1);
95 A(i, i - 2) += eim / h4;
96 A(i, i - 1) += -2.0 * (eim + eii) / h4;
97 A(i, i) += (eim + 4.0 * eii + eip) / h4 + beam_.k_at(i);
98 A(i, i + 1) += -2.0 * (eii + eip) / h4;
99 A(i, i + 2) += eip / h4;
100 b[i] = beam_.q_at(i);
101 }
102
103 // ---- Boundary rows -----------------------------------------------------
104 // The four stencils below are written for the left end in terms of the
105 // four boundary nodes (e0, e1, e2, e3) = (0, 1, 2, 3); the right end uses
106 // the mirror image (N, N-1, N-2, N-3) with the sign of the (odd-order)
107 // shear stencil flipped.
108 auto fill_end = [&](bool left_end) {
109 const Support sup = left_end ? beam_.left() : beam_.right();
110 // Node indices walking inward from the end, and the row indices for
111 // the outer (deflection/moment) and inner (slope/shear) conditions.
112 const std::size_t e0 = left_end ? 0 : N;
113 const std::size_t e1 = left_end ? 1 : N - 1;
114 const std::size_t e2 = left_end ? 2 : N - 2;
115 const std::size_t e3 = left_end ? 3 : N - 3;
116 const std::size_t outer_row = e0;
117 const std::size_t inner_row = e1;
118 const double ei0 = beam_.EI_at(e0);
119 const double ei1 = beam_.EI_at(e1);
120 const double ei2 = beam_.EI_at(e2);
121 // Forward (+1) at the left end, backward (-1) at the right end, used to
122 // orient the one-sided first-derivative (slope/shear) stencils.
123 const double dir = left_end ? 1.0 : -1.0;
124
125 // Second-derivative (w'') one-sided stencil at the end node:
126 // (2 w_e0 - 5 w_e1 + 4 w_e2 - w_e3) / h^2 (2nd order)
127 auto set_w2 = [&](std::size_t row, double scale, double rhs) {
128 A(row, e0) += scale * 2.0 / h2;
129 A(row, e1) += scale * -5.0 / h2;
130 A(row, e2) += scale * 4.0 / h2;
131 A(row, e3) += scale * -1.0 / h2;
132 b[row] = rhs;
133 };
134
135 // --- Outer row: deflection (essential) or moment (natural) ---------
136 if (sup == Support::Free) {
137 // EI * w'' = M_applied
138 set_w2(outer_row, ei0, beam_.end_moment(left_end));
139 } else {
140 // w = 0
141 A(outer_row, e0) += 1.0;
142 b[outer_row] = 0.0;
143 }
144
145 // --- Inner row: slope (clamped) / moment (pinned) / shear (free) ---
146 if (sup == Support::Clamped) {
147 // w' = 0, 2nd-order one-sided: (-3 w_e0 + 4 w_e1 - w_e2)/(2h) = 0
148 A(inner_row, e0) += dir * -3.0 / (2.0 * h);
149 A(inner_row, e1) += dir * 4.0 / (2.0 * h);
150 A(inner_row, e2) += dir * -1.0 / (2.0 * h);
151 b[inner_row] = 0.0;
152 } else if (sup == Support::Pinned) {
153 // EI * w'' = 0 (zero moment)
154 set_w2(inner_row, ei0, 0.0);
155 } else { // Free: (EI w'')' = V_applied
156 // Shear = M'(end) via a one-sided first difference of the nodal
157 // moment field M_j = EI_j w''_j:
158 // M_e0 = EI_e0 (2 w_e0 - 5 w_e1 + 4 w_e2 - w_e3)/h^2 (one-sided)
159 // M_e1 = EI_e1 ( w_e0 - 2 w_e1 + w_e2 )/h^2 (central)
160 // M_e2 = EI_e2 ( w_e1 - 2 w_e2 + w_e3 )/h^2 (central)
161 // M'(end) ~ (-3 M_e0 + 4 M_e1 - M_e2)/(2h)
162 // Combined coefficients (factor 1/(2 h^3)):
163 const double f = dir / (2.0 * h * h2);
164 A(inner_row, e0) += f * (-6.0 * ei0 + 4.0 * ei1);
165 A(inner_row, e1) += f * (15.0 * ei0 - 8.0 * ei1 - ei2);
166 A(inner_row, e2) += f * (-12.0 * ei0 + 4.0 * ei1 + 2.0 * ei2);
167 A(inner_row, e3) += f * (3.0 * ei0 - ei2);
168 b[inner_row] = beam_.applied_end_shear(left_end);
/opt/beamfd/src/solver.cpp
168 b[inner_row] = beam_.applied_end_shear(left_end);
169 }
170 };
171
172 fill_end(true);
173 fill_end(false);
174
175 return sys;
176 }
177
178 BeamResult BeamSolver::solve() const {
179 beam_.validate(); // propagate mechanism / ill-posed configuration errors
180
181 BeamSystem sys = assemble_system();
182 BeamResult result;
183 result.deflection = solve_general(sys.A, sys.b); // propagates if singular
184 return result;
185 }
186
187 // ----------------------------------------------------------------------------
188 // Recovery helpers (already implemented; do NOT modify). These differentiate a
189 // given deflection field so callers can recompute physical quantities (bending
190 // moment, shear) directly from a candidate result, honoring a variable EI(x).
191
192 namespace {
/opt/beamfd/src/solver.cpp
// ----------------------------------------------------------------------------
// Finite-difference discretization of
//
// d^2/dx^2 ( EI(x) w''(x) ) + k(x) w(x) = q(x)
//
// on the uniform grid x_i = i*dx, i = 0..N (n = N+1 nodes), second-order
// accurate throughout.
//
// The fourth-order operator is built by nesting two central second differences:
// the nodal moment is M_j = EI_j * w''_j with the *same* second-difference
// stencils used by the recovery helpers (central in the interior, a 2nd-order
// one-sided form at an end), and the governing equation is the second
// difference of that moment field,
//
// (M_{i-1} - 2 M_i + M_{i+1}) / dx^2 + k_i w_i = q_i,
//
// expanded into a five-point [w_{i-2} .. w_{i+2}] stencil with variable EI:
//
// w_{i-2}: EI_{i-1} \
// w_{i-1}: -2(EI_{i-1} + EI_i) | all /dx^4
// w_i : EI_{i-1} + 4 EI_i + EI_{i+1} |
// w_{i+1}: -2(EI_i + EI_{i+1}) |
// w_{i+2}: EI_{i+1} /
//
// This interior equation is enforced at nodes i = 2 .. N-2. The remaining four
// rows (0, 1, N-1, N) carry the two boundary conditions at each end:
//
// * "outer" row (0 / N): the deflection-or-moment condition
// clamped / pinned : w = 0
// free : EI*w'' = M_applied
// * "inner" row (1 / N-1): the slope-or-shear condition
// clamped : w' = 0
// pinned : EI*w'' = 0 (zero moment)
// free : (EI*w'')' = V_applied (shear)
//
// The slope, moment and shear stencils are all 2nd-order one-sided forms
// consistent with the interior scheme and with the recovery helpers.
BeamSystem BeamSolver::assemble_system() const {
const std::size_t n = beam_.num_nodes();
const std::size_t N = beam_.num_intervals(); // last node index
const double h = beam_.dx();
const double h2 = h * h;
const double h4 = h2 * h2;
BeamSystem sys;
sys.A = DenseMatrix(n, n);
sys.b.assign(n, 0.0);
DenseMatrix& A = sys.A;
std::vector<double>& b = sys.b;
// ---- Interior governing equation at nodes 2 .. N-2 ---------------------
for (std::size_t i = 2; i + 2 <= N; ++i) {
const double eim = beam_.EI_at(i - 1);
const double eii = beam_.EI_at(i);
const double eip = beam_.EI_at(i + 1);
A(i, i - 2) += eim / h4;
A(i, i - 1) += -2.0 * (eim + eii) / h4;
A(i, i) += (eim + 4.0 * eii + eip) / h4 + beam_.k_at(i);
A(i, i + 1) += -2.0 * (eii + eip) / h4;
A(i, i + 2) += eip / h4;
b[i] = beam_.q_at(i);
}
// ---- Boundary rows -----------------------------------------------------
// The four stencils below are written for the left end in terms of the
// four boundary nodes (e0, e1, e2, e3) = (0, 1, 2, 3); the right end uses
// the mirror image (N, N-1, N-2, N-3) with the sign of the (odd-order)
// shear stencil flipped.
auto fill_end = [&](bool left_end) {
const Support sup = left_end ? beam_.left() : beam_.right();
// Node indices walking inward from the end, and the row indices for
// the outer (deflection/moment) and inner (slope/shear) conditions.
const std::size_t e0 = left_end ? 0 : N;
const std::size_t e1 = left_end ? 1 : N - 1;
const std::size_t e2 = left_end ? 2 : N - 2;
const std::size_t e3 = left_end ? 3 : N - 3;
const std::size_t outer_row = e0;
const std::size_t inner_row = e1;
const double ei0 = beam_.EI_at(e0);
const double ei1 = beam_.EI_at(e1);
const double ei2 = beam_.EI_at(e2);
// Forward (+1) at the left end, backward (-1) at the right end, used to
// orient the one-sided first-derivative (slope/shear) stencils.
const double dir = left_end ? 1.0 : -1.0;
// Second-derivative (w'') one-sided stencil at the end node:
// (2 w_e0 - 5 w_e1 + 4 w_e2 - w_e3) / h^2 (2nd order)
auto set_w2 = [&](std::size_t row, double scale, double rhs) {
A(row, e0) += scale * 2.0 / h2;
A(row, e1) += scale * -5.0 / h2;
A(row, e2) += scale * 4.0 / h2;
A(row, e3) += scale * -1.0 / h2;
b[row] = rhs;
};
// --- Outer row: deflection (essential) or moment (natural) ---------
if (sup == Support::Free) {
// EI * w'' = M_applied
set_w2(outer_row, ei0, beam_.end_moment(left_end));
} else {
// w = 0
A(outer_row, e0) += 1.0;
b[outer_row] = 0.0;
}
// --- Inner row: slope (clamped) / moment (pinned) / shear (free) ---
if (sup == Support::Clamped) {
// w' = 0, 2nd-order one-sided: (-3 w_e0 + 4 w_e1 - w_e2)/(2h) = 0
A(inner_row, e0) += dir * -3.0 / (2.0 * h);
A(inner_row, e1) += dir * 4.0 / (2.0 * h);
A(inner_row, e2) += dir * -1.0 / (2.0 * h);
b[inner_row] = 0.0;
} else if (sup == Support::Pinned) {
// EI * w'' = 0 (zero moment)
set_w2(inner_row, ei0, 0.0);
} else { // Free: (EI w'')' = V_applied
// Shear = M'(end) via a one-sided first difference of the nodal
// moment field M_j = EI_j w''_j:
// M_e0 = EI_e0 (2 w_e0 - 5 w_e1 + 4 w_e2 - w_e3)/h^2 (one-sided)
// M_e1 = EI_e1 ( w_e0 - 2 w_e1 + w_e2 )/h^2 (central)
// M_e2 = EI_e2 ( w_e1 - 2 w_e2 + w_e3 )/h^2 (central)
// M'(end) ~ (-3 M_e0 + 4 M_e1 - M_e2)/(2h)
// Combined coefficients (factor 1/(2 h^3)):
const double f = dir / (2.0 * h * h2);
A(inner_row, e0) += f * (-6.0 * ei0 + 4.0 * ei1);
A(inner_row, e1) += f * (15.0 * ei0 - 8.0 * ei1 - ei2);
A(inner_row, e2) += f * (-12.0 * ei0 + 4.0 * ei1 + 2.0 * ei2);
A(inner_row, e3) += f * (3.0 * ei0 - ei2);
b[inner_row] = beam_.applied_end_shear(left_end);
}
};
fill_end(true);
fill_end(false);
return sys;
}// ----------------------------------------------------------------------------
// Finite-difference discretization of
//
// d^2/dx^2 ( EI(x) w''(x) ) + k(x) w(x) = q(x)
//
// on the uniform grid x_i = i*dx, i = 0..N (n = N+1 nodes), second-order
// accurate throughout (interior scheme, essential BCs and natural BCs alike).
//
// Nodal moment and the operator. Define the nodal moment with a central
// second difference, M_j = EI_j (w_{j-1} - 2 w_j + w_{j+1}) / dx^2, and take a
// second central difference of it for the governing operator,
//
// (M_{i-1} - 2 M_i + M_{i+1}) / dx^2 + k_i w_i = q_i,
//
// which expands into the variable-EI five-point [w_{i-2} .. w_{i+2}] stencil
//
// w_{i-2}: EI_{i-1} \
// w_{i-1}: -2(EI_{i-1} + EI_i) | all /dx^4
// w_i : EI_{i-1} + 4 EI_i + EI_{i+1} |
// w_{i+1}: -2(EI_i + EI_{i+1}) |
// w_{i+2}: EI_{i+1} /
//
// This governing equation is enforced at EVERY interior node i = 1 .. N-1
// (rows 1 .. N-1). At i = 1 it reaches the ghost node w_{-1} (through M_0) and
// at i = N-1 the ghost w_{N+1} (through M_N). Each ghost is expressed as a
// linear combination of nodal unknowns using ONE boundary condition per end,
//
// w_{-1} = c0 w_0 + c1 w_1 + g, (and the mirror at the right end)
//
// chosen so that the central second difference of the boundary nodal moment is
// consistent with that end's condition:
//
// clamped: w'(0)=0 (central) -> w_{-1} = w_1 (g = 0)
// pinned: w''(0)=0 -> w_{-1} = 2 w_0 - w_1 (g = 0)
// free: EI_0 w''(0)=M_app -> w_{-1} = 2 w_0 - w_1 + h^2 M_app/EI_0
//
// so that M_0 reproduces the prescribed end moment exactly (0 for clamped/
// pinned, M_app for free). The ghost is substituted into the i=1 equation, and
// its constant part g moves to the right-hand side.
//
// The two remaining rows (0 and N) carry the end's OTHER condition:
// clamped / pinned : w = 0 (essential)
// free : (EI w'')'(end) = V_app (shear, natural)
// The shear row uses a one-sided first difference of the nodal moment with the
// boundary moment fixed at its prescribed value M_app:
// left: (-3 M_0 + 4 M_1 - M_2)/(2h) = V_app, M_0 = M_app
// right: ( 3 M_N - 4 M_{N-1} + M_{N-2})/(2h) = V_app, M_N = M_app
//
// Keeping the governing equation at the near-boundary nodes (rather than
// overwriting them with a boundary stencil) is what preserves second-order
// accuracy for the natural (free-end) conditions.
BeamSystem BeamSolver::assemble_system() const {
const std::size_t n = beam_.num_nodes();
const std::size_t N = beam_.num_intervals(); // last node index
const double h = beam_.dx();
const double h2 = h * h;
const double h4 = h2 * h2;
BeamSystem sys;
sys.A = DenseMatrix(n, n);
sys.b.assign(n, 0.0);
DenseMatrix& A = sys.A;
std::vector<double>& b = sys.b;
// Ghost relation w_ghost = c_end * w_end + c_in * w_in + g for one end, with
// the boundary node "end" and its inward neighbour "in".
struct Ghost {
double c_end{0.0}; // coefficient on the boundary node
double c_in{0.0}; // coefficient on the inward neighbour
double g{0.0}; // constant term
};
auto ghost_for = [&](bool left_end) {
const Support sup = left_end ? beam_.left() : beam_.right();
const double ei_end = beam_.EI_at(left_end ? 0 : N);
Ghost gh;
switch (sup) {
case Support::Clamped: // w' = 0 (central): w_ghost = w_in
gh.c_end = 0.0;
gh.c_in = 1.0;
gh.g = 0.0;
break;
case Support::Pinned: // w'' = 0: w_ghost = 2 w_end - w_in
gh.c_end = 2.0;
gh.c_in = -1.0;
gh.g = 0.0;
break;
case Support::Free: // EI w'' = M_app: w_ghost = 2 w_end - w_in + h^2 M/EI
gh.c_end = 2.0;
gh.c_in = -1.0;
gh.g = h2 * beam_.end_moment(left_end) / ei_end;
break;
}
return gh;
};
const Ghost gl = ghost_for(true);
const Ghost gr = ghost_for(false);
// ---- Governing equation at interior nodes i = 1 .. N-1 ----------------
for (std::size_t i = 1; i + 1 <= N; ++i) {
const double eim = beam_.EI_at(i - 1); // EI at i-1 (left moment node)
const double eii = beam_.EI_at(i); // EI at i
const double eip = beam_.EI_at(i + 1); // EI at i+1 (right moment node)
// Standard five-point coefficients (some columns may be ghosts at the
// near-boundary nodes; they are folded in below).
double cm2 = eim / h4; // w_{i-2}
double cm1 = -2.0 * (eim + eii) / h4; // w_{i-1}
double c0 = (eim + 4.0 * eii + eip) / h4; // w_i
double cp1 = -2.0 * (eii + eip) / h4; // w_{i+1}
double cp2 = eip / h4; // w_{i+2}
b[i] = beam_.q_at(i) + beam_.k_at(i) * 0.0; // (k term added to c0 below)
c0 += beam_.k_at(i);
if (i == 1) {
// w_{i-2} == w_{-1} is the left ghost: substitute.
// contributes cm2 * (gl.c_end w_0 + gl.c_in w_1 + gl.g)
A(i, 0) += cm1 + cm2 * gl.c_end; // w_0 (== w_{i-1})
A(i, 1) += c0 + cm2 * gl.c_in; // w_1 (== w_i)
A(i, 2) += cp1; // w_2
A(i, 3) += cp2; // w_3
b[i] -= cm2 * gl.g;
} else if (i == N - 1) {
// w_{i+2} == w_{N+1} is the right ghost: substitute.
A(i, N) += cp1 + cp2 * gr.c_end; // w_N (== w_{i+1})
A(i, N - 1) += c0 + cp2 * gr.c_in; // w_{N-1} (== w_i)
A(i, N - 2) += cm1; // w_{N-2}
A(i, N - 3) += cm2; // w_{N-3}
b[i] -= cp2 * gr.g;
} else {
A(i, i - 2) += cm2;
A(i, i - 1) += cm1;
A(i, i) += c0;
A(i, i + 1) += cp1;
A(i, i + 2) += cp2;
}
}
// ---- Boundary rows (0 and N): essential w=0, or natural shear ---------
auto fill_boundary_row = [&](bool left_end) {
const Support sup = left_end ? beam_.left() : beam_.right();
const std::size_t e0 = left_end ? 0 : N; // boundary node
const std::size_t e1 = left_end ? 1 : N - 1; // first inward
const std::size_t e2 = left_end ? 2 : N - 2; // second inward
const std::size_t e3 = left_end ? 3 : N - 3; // third inward
if (sup != Support::Free) {
A(e0, e0) += 1.0; // w = 0
b[e0] = 0.0;
return;
}
// Free end: shear (EI w'')' = V_app, one-sided first difference of the
// nodal moment with the boundary moment held at the prescribed M_app:
// left : (-3 M_0 + 4 M_1 - M_2)/(2h) = V_app
// right: ( 3 M_N - 4 M_{N-1} + M_{N-2})/(2h) = V_app
// M_e0 = M_app (constant -> RHS); M_e1, M_e2 are central second diffs.
const double M_app = beam_.end_moment(left_end);
const double V_app = beam_.applied_end_shear(left_end);
const double ei1 = beam_.EI_at(e1);
const double ei2 = beam_.EI_at(e2);
const double s = 1.0 / (2.0 * h); // outer 1/(2h) of the shear difference
const double a1 = left_end ? 4.0 : -4.0; // coeff on M_{first inward}
const double a2 = left_end ? -1.0 : 1.0; // coeff on M_{second inward}
const double a0 = left_end ? -3.0 : 3.0; // coeff on M_end (-> RHS)
// a1 * M_e1 : EI_e1 (w_e0 - 2 w_e1 + w_e2)/h^2
A(e0, e0) += s * a1 * ei1 / h2;
A(e0, e1) += s * a1 * (-2.0) * ei1 / h2;
A(e0, e2) += s * a1 * ei1 / h2;
// a2 * M_e2 : EI_e2 (w_e1 - 2 w_e2 + w_e3)/h^2
A(e0, e1) += s * a2 * ei2 / h2;
A(e0, e2) += s * a2 * (-2.0) * ei2 / h2;
A(e0, e3) += s * a2 * ei2 / h2;
// a0 * M_e0 = a0 * M_app is constant -> move to RHS.
b[e0] = V_app - s * a0 * M_app;
};
fill_boundary_row(true);
fill_boundary_row(false);
return sys;
}/opt/beamfd/src/solver.cpp
double cp2 = eip / h4; // w_{i+2}
b[i] = beam_.q_at(i) + beam_.k_at(i) * 0.0; // (k term added to c0 below)
c0 += beam_.k_at(i); double cp2 = eip / h4; // w_{i+2}
c0 += beam_.k_at(i); // Winkler foundation term
b[i] = beam_.q_at(i);Rebuild and run all verification harnesses
cd /opt/beamfd && cmake --build build --parallel >/tmp/build.log 2>&1 && echo BUILD_OK && (cd build && ctest --output-on-failure 2>&1 | tail -4) echo "=== closed-form ==="; g++ -std=c++17 -Iinclude /tmp/verify.cpp src/beam.cpp src/solver.cpp src/linalg.cpp -o /tmp/verify && /tmp/verify echo "=== manufactured (varEI/foundation/free) ==="; g++ -std=c++17 -Iinclude /tmp/verify2.cpp src/beam.cpp src/solver.cpp src/linalg.cpp -o /tmp/verify2 && /tmp/verify2 echo "=== constEI free ==="; g++ -std=c++17 -Iinclude /tmp/verify3.cpp src/beam.cpp src/solver.cpp src/linalg.cpp -o /tmp/verify3 && /tmp/verify3
BUILD_OK 100% tests passed, 0 tests failed out of 1 Total Test time (real) = 0.00 sec === closed-form === pinned-pinned UDL N= 20 err=2.109e-03 N= 40 err=5.273e-04 rate=2.00 N= 80 err=1.318e-04 rate=2.00 N= 160 err=3.296e-05 rate=2.00 clamped-clamped UDL N= 20 err=4.219e-03 N= 40 err=1.055e-03 rate=2.00 N= 80 err=2.637e-04 rate=2.00 N= 160 err=6.592e-05 rate=2.00 cantilever tip load N= 20 err=1.125e-02 N= 40 err=2.813e-03 rate=2.00 N= 80 err=7.031e-04 rate=2.00 N= 160 err=1.759e-04 rate=2.00 === manufactured (varEI/foundation/free) === varEI clamped-clamped err=2.000e-02 err=5.000e-03(r=2.00) err=1.250e-03(r=2.00) err=3.125e-04(r=2.00) varEI+found clamp-clamp err=1.975e-02 err=4.939e-03(r=2.00) err=1.235e-03(r=2.00) err=3.087e-04(r=2.00) varEI clamped-FREE err=2.223e-01 err=5.546e-02(r=2.00) err=1.386e-02(r=2.00) err=3.464e-03(r=2.00) varEI+found clamped-FREE err=1.397e-01 err=3.478e-02(r=2.01) err=8.686e-03(r=2.00) err=2.171e-03(r=2.00) const EI+found c-c err=1.964e-02 err=4.911e-03(r=2.00) err=1.228e-03(r=2.00) err=3.069e-04(r=2.00) === constEI free === constEI clamped-FREE err=4.000e-02 err=1.000e-02(r=2.00) err=2.500e-03(r=2.00) err=6.250e-04(r=2.00) [stdout] BUILD_OK 100% tests passed, 0 tests failed out of 1 Total Test time (real) = 0.00 sec === closed-form === pinned-pinned UDL N= 20 err=2.109e-03 N= 40 err=5.273e-04 rate=2.00 N= 80 err=1.318e-04 rate=2.00 N= 160 err=3.296e-05 rate=2.00 clamped-clamped UDL N= 20 err=4.219e-03 N= 40 err=1.055e-03 rate=2.00 N= 80 err=2.637e-04 rate=2.00 N= 160 err=6.592e-05 rate=2.00 cantilever tip load N= 20 err=1.125e-02 N= 40 err=2.813e-03 rate=2.00 N= 80 err=7.031e-04 rate=2.00 N= 160 err=1.759e-04 rate=2.00 === manufactured (varEI/foundation/free) === varEI clamped-clamped err=2.000e-02 err=5.000e-03(r=2.00) err=1.250e-03(r=2.00) err=3.125e-04(r=2.00) varEI+found clamp-clamp err=1.975e-02 err=4.939e-03(r=2.00) err=1.235e-03(r=2.00) err=3.087e-04(r=2.00) varEI clamped-FREE err=2.223e-01 err=5.546e-02(r=2.00) err=1.386e-02(r=2.00) err=3.464e-03(r=2.00) varEI+found clamped-FREE err=1.397e-01 err=3.478e-02(r=2.01) err=8.686e-03(r=2.00) err=2.171e-03(r=2.00) const EI+found c-c err=1.964e-02 err=4.911e-03(r=2.00) err=1.228e-03(r=2.00) err=3.069e-04(r=2.00) === constEI free === constEI clamped-FREE err=4.000e-02 err=1.000e-02(r=2.00) err=2.500e-03(r=2.00) err=6.250e-04(r=2.00)
Test free-free and free-clamped with foundation
cat > /tmp/verify4.cpp <<'EOF'
#include <cmath>
#include <cstdio>
#include <vector>
#include "beamfd/beam.hpp"
#include "beamfd/solver.hpp"
using namespace beamfd;
// w = x^2 (L-x)^2 ; EI=p+qx ; k=c+dx . Compatible with Clamped or Free ends
// (we supply M,V at free ends). Always has a foundation if c or d > 0.
struct MS { double L,p,q,c,d;
double w(double x)const{return x*x*(L-x)*(L-x);}
double w2(double x)const{return 2*L*L-12*L*x+12*x*x;}
double EI(double x)const{return p+q*x;}
double k(double x)const{return c+d*x;}
double M(double x)const{return EI(x)*w2(x);}
double Mpp(double x)const{return 24*p-24*q*L+72*q*x;}
double Mp(double x)const{double g=w2(x),gp=-12*L+24*x;return q*g+EI(x)*gp;}
double load(double x)const{return Mpp(x)+k(x)*w(x);}
};
double run(const MS&ms, std::size_t N, Support L_, Support R_){
double L=ms.L; Beam b(L,1.0,N,L_,R_);
std::vector<double> ei(N+1),kk(N+1),q(N+1);
for(std::size_t i=0;i<=N;++i){double x=L*i/N; ei[i]=ms.EI(x); kk[i]=ms.k(x); q[i]=ms.load(x);}
b.set_ei_nodal(ei); b.set_foundation_nodal(kk); b.set_q_nodal(q);
if(L_==Support::Free){ b.set_end_moment(ms.M(0),true); b.set_end_load(-ms.Mp(0),true);}
if(R_==Support::Free){ b.set_end_moment(ms.M(L),false); b.set_end_load(-ms.Mp(L),false);}
auto w=BeamSolver(b).solve().deflection; double e=0;
for(std::size_t i=0;i<=N;++i){double x=L*i/N;e=std::max(e,std::fabs(w[i]-ms.w(x)));}
return e;
}
void conv(const char*name,const MS&ms,Support L_,Support R_){printf("%-26s ",name);double prev=0;
for(int k=0;k<4;++k){std::size_t N=20*(1<<k);double e=run(ms,N,L_,R_);printf("err=%.3e",e);
if(k)printf("(r=%.2f)",std::log(prev/e)/std::log(2.0));printf(" ");prev=e;}printf("\n");}
int main(){
MS f{2.0,1000.0,400.0,500.0,200.0}; // varEI + foundation (>0 so free-free ok)
conv("free-free +found", f, Support::Free, Support::Free);
conv("free-clamped +found",f, Support::Free, Support::Clamped);
MS g{2.0,1000.0,0.0,800.0,0.0}; // const EI + const foundation
conv("free-free constEI", g, Support::Free, Support::Free);
return 0;
}
EOF
g++ -std=c++17 -Iinclude /tmp/verify4.cpp src/beam.cpp src/solver.cpp src/linalg.cpp -o /tmp/verify4 && /tmp/verify4free-free +found err=4.624e-01 err=1.149e-01(r=2.01) err=2.869e-02(r=2.00) err=7.171e-03(r=2.00) free-clamped +found err=1.353e-01 err=3.382e-02(r=2.00) err=8.455e-03(r=2.00) err=2.114e-03(r=2.00) free-free constEI err=6.601e-03 err=1.635e-03(r=2.01) err=4.077e-04(r=2.00) err=1.018e-04(r=2.00) [stdout] free-free +found err=4.624e-01 err=1.149e-01(r=2.01) err=2.869e-02(r=2.00) err=7.171e-03(r=2.00) free-clamped +found err=1.353e-01 err=3.382e-02(r=2.00) err=8.455e-03(r=2.00) err=2.114e-03(r=2.00) free-free constEI err=6.601e-03 err=1.635e-03(r=2.01) err=4.077e-04(r=2.00) err=1.018e-04(r=2.00)
Test pinned-end variable-EI manufactured solutions
cat > /tmp/verify5.cpp <<'EOF'
#include <cmath>
#include <cstdio>
#include <vector>
#include "beamfd/beam.hpp"
#include "beamfd/solver.hpp"
using namespace beamfd;
// w=sin(a x), a=pi/L : w(0)=0,w''(0)=0 (pinned-compatible at x=0); supply M,V at free.
struct MS{ double L,p,q,c,d; double a()const{return M_PI/L;}
double w(double x)const{return std::sin(a()*x);}
double EI(double x)const{return p+q*x;}
double k(double x)const{return c+d*x;}
double M(double x)const{double A=a();return EI(x)*(-A*A*std::sin(A*x));}
double Mp(double x)const{double A=a();return q*(-A*A*std::sin(A*x))+EI(x)*(-A*A*A*std::cos(A*x));}
double Mpp(double x)const{double A=a();return -2*q*A*A*A*std::cos(A*x)+EI(x)*A*A*A*A*std::sin(A*x);}
double load(double x)const{return Mpp(x)+k(x)*w(x);}
};
double run(const MS&ms,std::size_t N,Support L_,Support R_){
double L=ms.L; Beam b(L,1.0,N,L_,R_);
std::vector<double> ei(N+1),kk(N+1),q(N+1);
for(std::size_t i=0;i<=N;++i){double x=L*i/N;ei[i]=ms.EI(x);kk[i]=ms.k(x);q[i]=ms.load(x);}
b.set_ei_nodal(ei);b.set_foundation_nodal(kk);b.set_q_nodal(q);
if(L_==Support::Free){b.set_end_moment(ms.M(0),true);b.set_end_load(-ms.Mp(0),true);}
if(R_==Support::Free){b.set_end_moment(ms.M(L),false);b.set_end_load(-ms.Mp(L),false);}
auto w=BeamSolver(b).solve().deflection;double e=0;
for(std::size_t i=0;i<=N;++i){double x=L*i/N;e=std::max(e,std::fabs(w[i]-ms.w(x)));}return e;}
void conv(const char*name,const MS&ms,Support L_,Support R_){printf("%-26s ",name);double prev=0;
for(int k=0;k<4;++k){std::size_t N=20*(1<<k);double e=run(ms,N,L_,R_);printf("err=%.3e",e);
if(k)printf("(r=%.2f)",std::log(prev/e)/std::log(2.0));printf(" ");prev=e;}printf("\n");}
int main(){
MS m{2.0,1500.0,600.0,0.0,0.0};
MS mf{2.0,1500.0,600.0,400.0,150.0};
conv("pinned-pinned varEI", m, Support::Pinned,Support::Pinned);
conv("pinned-pinned vEI+fnd",mf,Support::Pinned,Support::Pinned);
conv("pinned-free vEI+fnd", mf,Support::Pinned,Support::Free);
conv("free-pinned vEI+fnd", mf,Support::Free, Support::Pinned);
return 0;
}
EOF
g++ -std=c++17 -Iinclude /tmp/verify5.cpp src/beam.cpp src/solver.cpp src/linalg.cpp -o /tmp/verify5 && /tmp/verify5pinned-pinned varEI err=4.107e-03 err=1.025e-03(r=2.00) err=2.561e-04(r=2.00) err=6.404e-05(r=2.00) pinned-pinned vEI+fnd err=3.936e-03 err=9.823e-04(r=2.00) err=2.455e-04(r=2.00) err=6.138e-05(r=2.00) pinned-free vEI+fnd err=2.564e-01 err=6.589e-02(r=1.96) err=1.669e-02(r=1.98) err=4.201e-03(r=1.99) free-pinned vEI+fnd err=1.949e-01 err=4.658e-02(r=2.06) err=1.136e-02(r=2.04) err=2.806e-03(r=2.02) [stdout] pinned-pinned varEI err=4.107e-03 err=1.025e-03(r=2.00) err=2.561e-04(r=2.00) err=6.404e-05(r=2.00) pinned-pinned vEI+fnd err=3.936e-03 err=9.823e-04(r=2.00) err=2.455e-04(r=2.00) err=6.138e-05(r=2.00) pinned-free vEI+fnd err=2.564e-01 err=6.589e-02(r=1.96) err=1.669e-02(r=1.98) err=4.201e-03(r=1.99) free-pinned vEI+fnd err=1.949e-01 err=4.658e-02(r=2.06) err=1.136e-02(r=2.04) err=2.806e-03(r=2.02)
Test validate propagation, warnings, and demo CLI
cd /opt/beamfd
cat > /tmp/verify6.cpp <<'EOF'
#include <cstdio>
#include <stdexcept>
#include "beamfd/beam.hpp"
#include "beamfd/solver.hpp"
using namespace beamfd;
int main(){
// free-free without foundation -> validate() must throw out of solve()
try { Beam b(2.0,1000.0,20,Support::Free,Support::Free); b.set_distributed_load(1.0);
BeamSolver(b).solve(); printf("FAIL: no throw\n"); }
catch(const std::exception&e){ printf("OK validate threw: %s\n", e.what()); }
// free-pinned without foundation -> mechanism
try { Beam b(2.0,1000.0,20,Support::Free,Support::Pinned);
BeamSolver(b).solve(); printf("FAIL: no throw\n"); }
catch(const std::exception&e){ printf("OK validate threw: %s\n", e.what()); }
return 0;
}
EOF
g++ -std=c++17 -Wall -Wextra -Iinclude /tmp/verify6.cpp src/beam.cpp src/solver.cpp src/linalg.cpp -o /tmp/verify6 2>/tmp/warn.log && echo "--warnings--" && cat /tmp/warn.log && /tmp/verify6
echo "=== demo CLI on example ==="; ./build/* examples/cantilever_tip.beam 2>/dev/null | head -40 || ls build--warnings--
src/solver.cpp:55:1: warning: multi-line comment [-Wcomment]
55 | // w_{i-2}: EI_{i-1} \
| ^
OK validate threw: Beam::validate: free-free beam is a mechanism
OK validate threw: Beam::validate: a free end requires a clamped opposite end
=== demo CLI on example ===
[stdout]
--warnings--
src/solver.cpp:55:1: warning: multi-line comment [-Wcomment]
55 | // w_{i-2}: EI_{i-1} \
| ^
OK validate threw: Beam::validate: free-free beam is a mechanism
OK validate threw: Beam::validate: a free end requires a clamped opposite end
=== demo CLI on example ===/opt/beamfd/src/solver.cpp
// w_{i-2}: EI_{i-1} \
// w_{i-1}: -2(EI_{i-1} + EI_i) | all /dx^4
// w_i : EI_{i-1} + 4 EI_i + EI_{i+1} |
// w_{i+1}: -2(EI_i + EI_{i+1}) |
// w_{i+2}: EI_{i+1} /// w_{i-2}: EI_{i-1}
// w_{i-1}: -2(EI_{i-1} + EI_i) (all coefficients /dx^4)
// w_i : EI_{i-1} + 4 EI_i + EI_{i+1}
// w_{i+1}: -2(EI_i + EI_{i+1})
// w_{i+2}: EI_{i+1}Confirm clean build and find demo binary
g++ -std=c++17 -Wall -Wextra -Iinclude -c src/solver.cpp -o /tmp/s.o 2>&1 | head; echo "warnings above (none = clean)" ls build/ | head -20 echo "=== run demo ==="; find build -maxdepth 1 -type f -executable | head
warnings above (none = clean) CMakeCache.txt CMakeFiles CTestTestfile.cmake Makefile Testing beamfd_cli cmake_install.cmake libbeamfd.a tests === run demo === build/beamfd_cli [stdout] warnings above (none = clean) CMakeCache.txt CMakeFiles CTestTestfile.cmake Makefile Testing beamfd_cli cmake_install.cmake libbeamfd.a tests === run demo === build/beamfd_cli
Rebuild and run demo CLI plus ctest
cmake --build build --parallel >/tmp/b.log 2>&1 && echo REBUILT && ./build/beamfd_cli examples/cantilever_tip.beam 2>&1 | head -40 echo "=== ctest ==="; (cd build && ctest --output-on-failure 2>&1 | tail -4)
REBUILT # 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 -6.292651e-14 3.000000e+03 -1.000000e+03 1 1.500000e-02 1.687499e-07 2.985000e+03 -1.000000e+03 2 3.000000e-02 6.733124e-07 2.970000e+03 -1.000000e+03 3 4.500000e-02 1.512000e-06 2.955000e+03 -1.000000e+03 4 6.000000e-02 2.683125e-06 2.940000e+03 -1.000000e+03 5 7.500000e-02 4.185000e-06 2.925000e+03 -1.000000e+03 6 9.000000e-02 6.015937e-06 2.910000e+03 -1.000000e+03 7 1.050000e-01 8.174250e-06 2.895000e+03 -1.000000e+03 8 1.200000e-01 1.065825e-05 2.880000e+03 -1.000000e+03 9 1.350000e-01 1.346625e-05 2.865000e+03 -1.000000e+03 10 1.500000e-01 1.659656e-05 2.850000e+03 -1.000000e+03 11 1.650000e-01 2.004750e-05 2.835000e+03 -1.000000e+03 12 1.800000e-01 2.381737e-05 2.820000e+03 -1.000000e+03 13 1.950000e-01 2.790450e-05 2.805000e+03 -1.000000e+03 14 2.100000e-01 3.230719e-05 2.790000e+03 -1.000000e+03 15 2.250000e-01 3.702375e-05 2.775000e+03 -1.000000e+03 16 2.400000e-01 4.205250e-05 2.760000e+03 -1.000000e+03 17 2.550000e-01 4.739175e-05 2.745000e+03 -1.000000e+03 18 2.700000e-01 5.303981e-05 2.730000e+03 -1.000000e+03 19 2.850000e-01 5.899500e-05 2.715000e+03 -1.000000e+03 20 3.000000e-01 6.525562e-05 2.700000e+03 -1.000000e+03 21 3.150000e-01 7.182000e-05 2.685000e+03 -1.000000e+03 22 3.300000e-01 7.868644e-05 2.670000e+03 -1.000000e+03 23 3.450000e-01 8.585325e-05 2.655000e+03 -1.000000e+03 24 3.600000e-01 9.331875e-05 2.640000e+03 -1.000000e+03 25 3.750000e-01 1.010812e-04 2.625000e+03 -1.000000e+03 26 3.900000e-01 1.091391e-04 2.610000e+03 -1.000000e+03 27 4.050000e-01 1.174905e-04 2.595000e+03 -1.000000e+03 28 4.200000e-01 1.261339e-04 2.580000e+03 -1.000000e+03 29 4.350000e-01 1.350675e-04 2.565000e+03 -1.000000e+03 30 4.500000e-01 1.442897e-04 2.550000e+03 -1.000000e+03 31 4.650000e-01 1.537987e-04 2.535000e+03 -1.000000e+03 32 4.800000e-01 1.635930e-04 2.520000e+03 -1.000000e+03 33 4.950000e-01 1.736707e-04 2.505000e+03 -1.000000e+03 34 5.100000e-01 1.840303e-04 2.490000e+03 -1.000000e+03 35 5.250000e-01 1.946700e-04 2.475000e+03 -1.000000e+03 === ctest === 100% tests passed, 0 tests failed out of 1 Total Test time (real) = 0.00 sec [stdout] REBUILT # 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 -6.292651e-14 3.000000e+03 -1.000000e+03 1 1.500000e-02 1.687499e-07 2.985000e+03 -1.000000e+03 2 3.000000e-02 6.733124e-07 2.970000e+03 -1.000000e+03 3 4.500000e-02 1.512000e-06 2.955000e+03 -1.000000e+03 4 6.000000e-02 2.683125e-06 2.940000e+03 -1.000000e+03 5 7.500000e-02 4.185000e-06 2.925000e+03 -1.000000e+03 6 9.000000e-02 6.015937e-06 2.910000e+03 -1.000000e+03 7 1.050000e-01 8.174250e-06 2.895000e+03 -1.000000e+03 8 1.200000e-01 1.065825e-05 2.880000e+03 -1.000000e+03 9 1.350000e-01 1.346625e-05 2.865000e+03 -1.000000e+03 10 1.500000e-01 1.659656e-05 2.850000e+03 -1.000000e+03 11 1.650000e-01 2.004750e-05 2.835000e+03 -1.000000e+03 12 1.800000e-01 2.381737e-05 2.820000e+03 -1.000000e+03 13 1.950000e-01 2.790450e-05 2.805000e+03 -1.000000e+03 14 2.100000e-01 3.230719e-05 2.790000e+03 -1.000000e+03 15 2.250000e-01 3.702375e-05 2.775000e+03 -1.000000e+03 16 2.400000e-01 4.205250e-05 2.760000e+03 -1.000000e+03 17 2.550000e-01 4.739175e-05 2.745000e+03 -1.000000e+03 18 2.700000e-01 5.303981e-05 2.730000e+03 -1.000000e+03 19 2.850000e-01 5.899500e-05 2.715000e+03 -1.000000e+03 20 3.000000e-01 6.525562e-05 2.700000e+03 -1.000000e+03 21 3.150000e-01 7.182000e-05 2.685000e+03 -1.000000e+03 22 3.300000e-01 7.868644e-05 2.670000e+03 -1.000000e+03 23 3.450000e-01 8.585325e-05 2.655000e+03 -1.000000e+03 24 3.600000e-01 9.331875e-05 2.640000e+03 -1.000000e+03 25 3.750000e-01 1.010812e-04 2.625000e+03 -1.000000e+03 26 3.900000e-01 1.091391e-04 2.610000e+03 -1.000000e+03 27 4.050000e-01 1.174905e-04 2.595000e+03 -1.000000e+03 28 4.200000e-01 1.261339e-04 2.580000e+03 -1.000000e+03 29 4.350000e-01 1.350675e-04 2.565000e+03 -1.000000e+03 30 4.500000e-01 1.442897e-04 2.550000e+03 -1.000000e+03 31 4.650000e-01 1.537987e-04 2.535000e+03 -1.000000e+03 32 4.800000e-01 1.635930e-04 2.520000e+03 -1.000000e+03 33 4.950000e-01 1.736707e-04 2.505000e+03 -1.000000e+03 34 5.100000e-01 1.840303e-04 2.490000e+03 -1.000000e+03 35 5.250000e-01 1.946700e-04 2.475000e+03 -1.000000e+03 === ctest === 100% tests passed, 0 tests failed out of 1 Total Test time (real) = 0.00 sec
/opt/beamfd/src/solver.cpp
108 double c_in{0.0}; // coefficient on the inward neighbour
109 double g{0.0}; // constant term
110 };
111 auto ghost_for = [&](bool left_end) {
112 const Support sup = left_end ? beam_.left() : beam_.right();
113 const double ei_end = beam_.EI_at(left_end ? 0 : N);
114 Ghost gh;
115 switch (sup) {
116 case Support::Clamped: // w' = 0 (central): w_ghost = w_in
117 gh.c_end = 0.0;
118 gh.c_in = 1.0;
119 gh.g = 0.0;
120 break;
121 case Support::Pinned: // w'' = 0: w_ghost = 2 w_end - w_in
122 gh.c_end = 2.0;
123 gh.c_in = -1.0;
124 gh.g = 0.0;
125 break;
126 case Support::Free: // EI w'' = M_app: w_ghost = 2 w_end - w_in + h^2 M/EI
127 gh.c_end = 2.0;
128 gh.c_in = -1.0;
129 gh.g = h2 * beam_.end_moment(left_end) / ei_end;
130 break;
131 }
132 return gh;
133 };
134
135 const Ghost gl = ghost_for(true);
136 const Ghost gr = ghost_for(false);
137
138 // ---- Governing equation at interior nodes i = 1 .. N-1 ----------------
139 for (std::size_t i = 1; i + 1 <= N; ++i) {
140 const double eim = beam_.EI_at(i - 1); // EI at i-1 (left moment node)
141 const double eii = beam_.EI_at(i); // EI at i
142 const double eip = beam_.EI_at(i + 1); // EI at i+1 (right moment node)
143
144 // Standard five-point coefficients (some columns may be ghosts at the
145 // near-boundary nodes; they are folded in below).
146 double cm2 = eim / h4; // w_{i-2}
147 double cm1 = -2.0 * (eim + eii) / h4; // w_{i-1}
148 double c0 = (eim + 4.0 * eii + eip) / h4; // w_i
149 double cp1 = -2.0 * (eii + eip) / h4; // w_{i+1}
150 double cp2 = eip / h4; // w_{i+2}
151 c0 += beam_.k_at(i); // Winkler foundation term
152 b[i] = beam_.q_at(i);
153
154 if (i == 1) {
155 // w_{i-2} == w_{-1} is the left ghost: substitute.
156 // contributes cm2 * (gl.c_end w_0 + gl.c_in w_1 + gl.g)
157 A(i, 0) += cm1 + cm2 * gl.c_end; // w_0 (== w_{i-1})
158 A(i, 1) += c0 + cm2 * gl.c_in; // w_1 (== w_i)
159 A(i, 2) += cp1; // w_2
160 A(i, 3) += cp2; // w_3
161 b[i] -= cm2 * gl.g;
162 } else if (i == N - 1) {
163 // w_{i+2} == w_{N+1} is the right ghost: substitute.
164 A(i, N) += cp1 + cp2 * gr.c_end; // w_N (== w_{i+1})
165 A(i, N - 1) += c0 + cp2 * gr.c_in; // w_{N-1} (== w_i)
166 A(i, N - 2) += cm1; // w_{N-2}
167 A(i, N - 3) += cm2; // w_{N-3}
168 b[i] -= cp2 * gr.g;
169 } else {
170 A(i, i - 2) += cm2;
171 A(i, i - 1) += cm1;
172 A(i, i) += c0;
173 A(i, i + 1) += cp1;
174 A(i, i + 2) += cp2;
175 }
176 }
177
178 // ---- Boundary rows (0 and N): essential w=0, or natural shear ---------
179 auto fill_boundary_row = [&](bool left_end) {
180 const Support sup = left_end ? beam_.left() : beam_.right();
181 const std::size_t e0 = left_end ? 0 : N; // boundary node
182 const std::size_t e1 = left_end ? 1 : N - 1; // first inward
183 const std::size_t e2 = left_end ? 2 : N - 2; // second inward
184 const std::size_t e3 = left_end ? 3 : N - 3; // third inward
185
186 if (sup != Support::Free) {
187 A(e0, e0) += 1.0; // w = 0
188 b[e0] = 0.0;
189 return;
190 }
191
192 // Free end: shear (EI w'')' = V_app, one-sided first difference of the
193 // nodal moment with the boundary moment held at the prescribed M_app:
194 // left : (-3 M_0 + 4 M_1 - M_2)/(2h) = V_app
195 // right: ( 3 M_N - 4 M_{N-1} + M_{N-2})/(2h) = V_app
196 // M_e0 = M_app (constant -> RHS); M_e1, M_e2 are central second diffs.
197 const double M_app = beam_.end_moment(left_end);
198 const double V_app = beam_.applied_end_shear(left_end);
199 const double ei1 = beam_.EI_at(e1);
200 const double ei2 = beam_.EI_at(e2);
201 const double s = 1.0 / (2.0 * h); // outer 1/(2h) of the shear difference
202 const double a1 = left_end ? 4.0 : -4.0; // coeff on M_{first inward}
Clean rebuild and full ctest
rm -rf build && cmake -S . -B build -DCMAKE_BUILD_TYPE=Release >/dev/null 2>&1 && cmake --build build --parallel >/tmp/b.log 2>&1 && echo "BUILD OK" && (cd build && ctest --output-on-failure 2>&1 | tail -5)
BUILD OK 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 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
No file changes, the agent left the workspace untouched.
== Configuring (cmake) ==
-- The CXX compiler identification is GNU 11.4.0
-- Detecting CXX compiler ABI info
-- Detecting CXX compiler ABI info - done
-- Check for working CXX compiler: /usr/bin/c++ - skipped
-- Detecting CXX compile features
-- Detecting CXX compile features - done
-- Configuring done
-- Generating done
-- Build files have been written to: /tmp/tmp.FW4FQxWH8k/beamfd/build_grade
== Building ==
[ 6%] Building CXX object CMakeFiles/beamfd.dir/src/linalg.cpp.o
[ 13%] Building CXX object CMakeFiles/beamfd.dir/src/solver.cpp.o
[ 20%] Building CXX object CMakeFiles/beamfd.dir/src/beam.cpp.o
[ 26%] Building CXX object CMakeFiles/beamfd.dir/src/io.cpp.o
[ 33%] Linking CXX static library libbeamfd.a
[ 33%] Built target beamfd
[ 40%] Building CXX object CMakeFiles/beamfd_cli.dir/src/main.cpp.o
[ 46%] Building CXX object tests/CMakeFiles/test_mms.dir/test_mms.cpp.o
[ 53%] Building CXX object tests/CMakeFiles/test_invariants.dir/test_invariants.cpp.o
[ 60%] Building CXX object tests/CMakeFiles/test_soak.dir/test_soak.cpp.o
[ 66%] Building CXX object tests/CMakeFiles/test_convergence.dir/test_convergence.cpp.o
[ 73%] Linking CXX executable beamfd_cli
[ 73%] Built target beamfd_cli
[ 80%] Linking CXX executable test_mms
[ 86%] Linking CXX executable test_convergence
[ 86%] Built target test_mms
[ 93%] Linking CXX executable test_invariants
[100%] Linking CXX executable test_soak
[100%] Built target test_convergence
[100%] Built target test_invariants
[100%] Built target test_soak
== Running hidden test suite (ctest) ==
Test project /tmp/tmp.FW4FQxWH8k/beamfd/build_grade
Start 1: test_mms
1/4 Test #1: test_mms ......................... Passed 0.01 sec
Start 2: test_convergence
2/4 Test #2: test_convergence ................. Passed 0.00 sec
Start 3: test_invariants
3/4 Test #3: test_invariants .................. Passed 0.00 sec
Start 4: test_soak
4/4 Test #4: test_soak ........................ Passed 0.01 sec
100% tests passed, 0 tests failed out of 4
Label Time Summary:
hidden = 0.02 sec*proc (4 tests)
Total Test time (real) = 0.03 sec
PASS: all hidden tests passedReproduce this trial: git checkout 2f94510 && PYTHONPATH=src python3 scripts/build_site.py , then open trial/trial_32b9f54b7dac41eb. Re-running the agent live requires EVAL_PLATFORM_ENABLE_OAUTH_SMOKE=1 and is non-deterministic.
Trial trial_32b9f54b7dac41eb · verifier authoritative; classifier explanatory.