tests/test.sh exit code (0 → resolved); the classification below is post-hoc and cannot change it./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/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/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
find . -type f -not -path './build/*' -not -path './.git/*' | sort
./CMakeLists.txt ./examples/cantilever_tip.beam ./include/beamfd/beam.hpp ./include/beamfd/io.hpp ./include/beamfd/linalg.hpp ./include/beamfd/solver.hpp ./include/beamfd/types.hpp ./src/beam.cpp ./src/io.cpp ./src/linalg.cpp ./src/main.cpp ./src/solver.cpp ./tests/CMakeLists.txt ./tests/test_smoke.cpp ./tests/test_util.hpp [stdout] ./CMakeLists.txt ./examples/cantilever_tip.beam ./include/beamfd/beam.hpp ./include/beamfd/io.hpp ./include/beamfd/linalg.hpp ./include/beamfd/solver.hpp ./include/beamfd/types.hpp ./src/beam.cpp ./src/io.cpp ./src/linalg.cpp ./src/main.cpp ./src/solver.cpp ./tests/CMakeLists.txt ./tests/test_smoke.cpp ./tests/test_util.hpp
/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
/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/src/io.cpp
1 #include "beamfd/io.hpp"
2
3 #include <array>
4 #include <fstream>
5 #include <iomanip>
6 #include <optional>
7 #include <sstream>
8 #include <stdexcept>
9 #include <string>
10 #include <utility>
11 #include <vector>
12
13 namespace beamfd {
14
15 namespace {
16
17 [[noreturn]] void fail(std::size_t line_no, const std::string& msg) {
18 std::ostringstream oss;
19 oss << "parse_beam: line " << line_no << ": " << msg;
20 throw std::runtime_error(oss.str());
21 }
22
23 Support parse_support(std::size_t line_no, const std::string& tok) {
24 if (tok == "clamped") return Support::Clamped;
25 if (tok == "pinned") return Support::Pinned;
26 if (tok == "free") return Support::Free;
27 fail(line_no, "support must be one of: clamped | pinned | free (got '" + tok +
28 "')");
29 }
30
31 bool parse_side(std::size_t line_no, const std::string& tok) {
32 if (tok == "left") return true;
33 if (tok == "right") return false;
34 fail(line_no, "side must be 'left' or 'right'");
35 }
36
37 std::vector<std::pair<double, double>> read_pairs(std::size_t line_no,
38 std::istringstream& ls,
39 const char* what) {
40 std::vector<std::pair<double, double>> pts;
41 double x, v;
42 while (ls >> x >> v) pts.emplace_back(x, v);
43 if (pts.empty()) fail(line_no, std::string("expected x value pairs after ") + what);
44 return pts;
45 }
46
47 } // namespace
48
49 Beam parse_beam(std::istream& in) {
50 std::optional<double> length, EI;
51 std::optional<std::size_t> cells;
52 std::optional<Support> left, right;
53 std::vector<std::pair<double, double>> ei_profile, foundation_profile;
54 std::optional<double> foundation_const;
55 double q = 0.0;
56 bool have_udl = false;
57 std::vector<std::array<double, 4>> loads;
58 std::vector<std::pair<bool, double>> end_forces, end_moments;
59
60 std::string line;
61 std::size_t line_no = 0;
62 while (std::getline(in, line)) {
63 ++line_no;
64 const auto hash = line.find('#');
65 if (hash != std::string::npos) line.erase(hash);
66 std::istringstream ls(line);
67 std::string tag;
68 if (!(ls >> tag)) continue;
69
70 if (tag == "length") {
71 double v;
72 if (!(ls >> v)) fail(line_no, "expected: length <L>");
73 length = v;
74 } else if (tag == "ei") {
75 double v;
76 if (!(ls >> v)) fail(line_no, "expected: ei <EI>");
77 EI = v;
78 } else if (tag == "ei_profile") {
79 ei_profile = read_pairs(line_no, ls, "ei_profile");
80 } else if (tag == "cells") {
81 std::size_t v;
82 if (!(ls >> v)) fail(line_no, "expected: cells <n>");
83 cells = v;
84 } else if (tag == "support") {
85 std::string a, b;
86 if (!(ls >> a >> b)) fail(line_no, "expected: support <left> <right>");
87 left = parse_support(line_no, a);
88 right = parse_support(line_no, b);
89 } else if (tag == "foundation") {
90 double v;
91 if (!(ls >> v)) fail(line_no, "expected: foundation <k>");
92 foundation_const = v;
93 } else if (tag == "foundation_profile") {
94 foundation_profile = read_pairs(line_no, ls, "foundation_profile");
95 } else if (tag == "udl") {
96 double v;
97 if (!(ls >> v)) fail(line_no, "expected: udl <q>");
98 q = v;
99 have_udl = true;
100 } else if (tag == "load") {
101 double x0, x1, q0, q1;
102 if (!(ls >> x0 >> x1 >> q0 >> q1))
103 fail(line_no, "expected: load <x0> <x1> <q0> <q1>");
104 loads.push_back({x0, x1, q0, q1});
105 } else if (tag == "endload") {
106 std::string side;
107 double v;
108 if (!(ls >> side >> v)) fail(line_no, "expected: endload <left|right> <P>");
109 end_forces.emplace_back(parse_side(line_no, side), v);
110 } else if (tag == "endmoment") {
111 std::string side;
112 double v;
113 if (!(ls >> side >> v))
114 fail(line_no, "expected: endmoment <left|right> <M>");
115 end_moments.emplace_back(parse_side(line_no, side), v);
116 } else {
117 fail(line_no, "unknown record '" + tag + "'");
118 }
119 }
120
121 if (!length) fail(line_no, "missing required record: length");
122 if (!cells) fail(line_no, "missing required record: cells");
123 if (!left || !right) fail(line_no, "missing required record: support");
124 if (!EI && ei_profile.empty())
125 fail(line_no, "missing required record: ei or ei_profile");
126
127 const double ei_seed = EI ? *EI : 1.0;
128 Beam beam(*length, ei_seed, *cells, *left, *right);
129 if (!ei_profile.empty()) beam.set_ei_profile(ei_profile);
130 if (foundation_const) beam.set_foundation_nodal(
131 std::vector<double>(beam.num_nodes(), *foundation_const));
132 if (!foundation_profile.empty()) beam.set_foundation_profile(foundation_profile);
133 if (have_udl) beam.set_distributed_load(q);
134 for (const auto& s : loads) beam.add_load_segment(s[0], s[1], s[2], s[3]);
135
136 try {
137 for (const auto& e : end_forces) beam.set_end_load(e.second, e.first);
138 for (const auto& e : end_moments) beam.set_end_moment(e.second, e.first);
139 } catch (const std::exception& ex) {
140 fail(line_no, ex.what());
141 }
142 return beam;
143 }
144
145 Beam parse_beam_file(const std::string& path) {
146 std::ifstream in(path);
147 if (!in) {
148 throw std::runtime_error("parse_beam_file: cannot open '" + path + "'");
149 }
150 return parse_beam(in);
151 }
152
153 void write_report(std::ostream& out, const Beam& beam, const BeamResult& result) {
154 BeamSolver solver(beam);
155 const std::vector<double> M = solver.bending_moment_field(result.deflection);
156 const std::vector<double> V = solver.shear_field(result.deflection);
157
158 out << std::scientific << std::setprecision(6);
159 out << "# beam analysis report\n";
160 out << "length=" << beam.length() << " EI(max)=" << beam.EI()
161 << " nodes=" << beam.num_nodes() << " dx=" << beam.dx() << "\n\n";
162
163 out << "[fields] (i, x, w[m], M=EI*w''[N*m], V=(EI*w'')'[N])\n";
164 for (std::size_t i = 0; i < beam.num_nodes(); ++i) {
165 out << i << " " << beam.node_x(i) << " " << result.deflection[i] << " "
166 << M[i] << " " << V[i] << "\n";
167 }
168 }
169
170 } // namespace beamfd
171
/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
//
// (EI(x) w''(x))'' + k(x) w(x) = q(x)
//
// on the uniform grid x_i = i*dx, i = 0..N, with n = N+1 nodal unknowns.
//
// Interior nodes (i = 2 .. N-2) use the conservative ("divergence-form") second-
// order 5-point stencil obtained by nesting two central second differences:
//
// M_i = EI_i (w_{i-1} - 2 w_i + w_{i+1}) / h^2 (curvature)
// M''_i = (M_{i-1} - 2 M_i + M_{i+1}) / h^2
//
// which is exact for cubics and O(h^2) in general. The four remaining equations
// (two per end) impose the boundary conditions to the same second order, using
// exactly the one-sided difference stencils that the recovery helpers above use
// for M = EI w'' and V = (EI w'')', so a recovered field reproduces the imposed
// moment/shear. Essential conditions (w, w') are imposed directly; natural ones
// (M, V at a free end) match the applied end actions.
// ----------------------------------------------------------------------------
BeamSystem BeamSolver::assemble_system() const {
const std::size_t n = beam_.num_nodes();
if (n < 5) {
// The second-order one-sided boundary stencils (and the shear recovery)
// need at least four nodes adjacent to each end; the framework's
// shear_field carries the same minimum.
throw std::runtime_error("BeamSolver::assemble_system: need at least 5 nodes");
}
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;
// Add coef * EI(j) * w''(x_j) (using the same second-difference stencil as
// second_derivative: one-sided at the ends, central in the interior) to the
// matrix row r. The 1/h^2 of the second difference is folded in here.
auto add_curvature = [&](std::size_t r, std::size_t j, double coef) {
const double f = coef * beam_.EI_at(j) / h2;
if (j == 0) {
A(r, 0) += 2.0 * f;
A(r, 1) += -5.0 * f;
A(r, 2) += 4.0 * f;
A(r, 3) += -1.0 * f;
} else if (j == n - 1) {
A(r, n - 1) += 2.0 * f;
A(r, n - 2) += -5.0 * f;
A(r, n - 3) += 4.0 * f;
A(r, n - 4) += -1.0 * f;
} else {
A(r, j - 1) += 1.0 * f;
A(r, j) += -2.0 * f;
A(r, j + 1) += 1.0 * f;
}
};
// ---- Interior equilibrium: (EI w'')'' + k w = q (rows i = 2 .. N-2) ------
for (std::size_t i = 2; i + 2 < n; ++i) {
const double eim = beam_.EI_at(i - 1);
const double ei = 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 - 2.0 * ei) / h4;
A(i, i) += (eim + 4.0 * ei + eip) / h4 + beam_.k_at(i);
A(i, i + 1) += (-2.0 * ei - 2.0 * eip) / h4;
A(i, i + 2) += eip / h4;
sys.b[i] = beam_.q_at(i);
}
// ---- Boundary rows -------------------------------------------------------
// Each end contributes its two conditions. The left end fills rows {0, 1}
// (which reference low-index nodes), the right end fills rows {n-1, n-2}.
// w = 0 at node e.
auto set_deflection_zero = [&](std::size_t r, std::size_t e) {
A(r, e) += 1.0;
sys.b[r] = 0.0;
};
// w' = 0 at an end, via a second-order one-sided first difference.
auto set_slope_zero = [&](std::size_t r, bool left) {
if (left) {
A(r, 0) += -3.0;
A(r, 1) += 4.0;
A(r, 2) += -1.0;
} else {
A(r, n - 1) += 3.0;
A(r, n - 2) += -4.0;
A(r, n - 3) += 1.0;
}
sys.b[r] = 0.0;
};
// EI w'' = M at node e (M = 0 for a pinned support).
auto set_moment = [&](std::size_t r, std::size_t e, double M) {
add_curvature(r, e, 1.0);
sys.b[r] = M;
};
// (EI w'')' = V at an end, via a one-sided first difference of the nodal
// moment field, matching shear_field.
auto set_shear = [&](std::size_t r, bool left, double V) {
const double inv2h = 1.0 / (2.0 * h);
if (left) {
add_curvature(r, 0, -3.0 * inv2h);
add_curvature(r, 1, 4.0 * inv2h);
add_curvature(r, 2, -1.0 * inv2h);
} else {
add_curvature(r, n - 1, 3.0 * inv2h);
add_curvature(r, n - 2, -4.0 * inv2h);
add_curvature(r, n - 3, 1.0 * inv2h);
}
sys.b[r] = V;
};
auto apply_end = [&](bool left) {
const Support s = left ? beam_.left() : beam_.right();
const std::size_t e = left ? 0u : n - 1;
const std::size_t r_first = left ? 0u : n - 1; // "primary" condition
const std::size_t r_second = left ? 1u : n - 2; // "secondary" condition
switch (s) {
case Support::Clamped:
set_deflection_zero(r_first, e);
set_slope_zero(r_second, left);
break;
case Support::Pinned:
set_deflection_zero(r_first, e);
set_moment(r_second, e, 0.0);
break;
case Support::Free:
set_moment(r_first, e, beam_.end_moment(left));
set_shear(r_second, left, beam_.applied_end_shear(left));
break;
}
};
apply_end(true);
apply_end(false);
return sys;
}
BeamResult BeamSolver::solve() const {
beam_.validate(); // propagate mechanism / ill-posed configuration errors
const BeamSystem sys = assemble_system();
BeamResult result;
// Propagates std::runtime_error from solve_general if A is singular.
result.deflection = solve_general(sys.A, sys.b);
return result;
}cmake -S . -B build -DCMAKE_BUILD_TYPE=Release >/tmp/cm.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 sec#include <cmath>
#include <cstdio>
#include <vector>
#include "beamfd/beam.hpp"
#include "beamfd/solver.hpp"
using namespace beamfd;
static double maxerr(const std::vector<double>&a,const std::vector<double>&b){
double m=0; for(size_t i=0;i<a.size();++i) m=std::max(m,std::fabs(a[i]-b[i])); return m;}
// Case A: cubic exact, clamped-free cantilever, tip shear only.
void caseA(){
double L=2.0, EI=1000.0; size_t N=40;
Beam beam(L,EI,N,Support::Clamped,Support::Free);
// manufactured w = 3L x^2 - x^3 -> w(0)=0,w'(0)=0; M(L)=0; V=(EI w'')'=EI*(-6)
// applied_end_shear(false) must equal -6EI -> -P = -6EI -> P=6EI
beam.set_end_load(6.0*EI,false);
beam.set_end_moment(0.0,false);
BeamSolver s(beam); auto r=s.solve();
std::vector<double> ex(N+1);
for(size_t i=0;i<=N;++i){double x=beam.node_x(i); ex[i]=3*L*x*x-x*x*x;}
printf("CaseA cantilever cubic: maxerr=%.3e (w(L)=%.6f exp %.6f)\n",
maxerr(r.deflection,ex), r.deflection[N], ex[N]);
}
// Case B: pinned-pinned UDL constant EI, w_max=5qL^4/384EI
void caseB(){
double L=3.0,EI=500.0,q=10.0; size_t N=80;
Beam beam(L,EI,N,Support::Pinned,Support::Pinned);
beam.set_distributed_load(q);
BeamSolver s(beam); auto r=s.solve();
double wmax=0; for(double w:r.deflection) wmax=std::max(wmax,std::fabs(w));
double exact=5*q*L*L*L*L/(384*EI);
printf("CaseB pinned UDL: wmax=%.6e exact=%.6e rel=%.3e\n",wmax,exact,std::fabs(wmax-exact)/exact);
}
// Case C: clamped-clamped, variable EI, foundation, manufactured w=x^2(L-x)^2
void caseC(size_t N,double&err){
double L=2.0,E0=800.0,b=0.5/L,k0=300.0;
Beam beam(L,E0,N,Support::Clamped,Support::Clamped);
std::vector<double> ei(N+1),kk(N+1),qq(N+1),ex(N+1);
for(size_t i=0;i<=N;++i){
double x=beam.node_x(i);
double w=x*x*(L-x)*(L-x);
double w2=2*L*L-12*L*x+12*x*x;
double w3=-12*L+24*x;
double EIx=E0*(1+b*x);
double M2=2*E0*b*w3 + 24*E0*(1+b*x);
ei[i]=EIx; kk[i]=k0; qq[i]=M2+k0*w; ex[i]=w;
}
beam.set_ei_nodal(ei); beam.set_foundation_nodal(kk); beam.set_q_nodal(qq);
BeamSolver s(beam); auto r=s.solve();
err=maxerr(r.deflection,ex);
}
// Case D: convergence for free-free on foundation (manufactured), variable EI
void caseD(size_t N,double&err){
double L=2.0,E0=600.0,b=0.3/L,k0=1000.0;
Beam beam(L,E0,N,Support::Free,Support::Free);
std::vector<double> ei(N+1),kk(N+1),qq(N+1),ex(N+1);
// w = cos(pi x/L) (nonzero at ends, slopes zero at ends but moment/shear nonzero)
double pi=std::acos(-1.0), w0=pi/L;
for(size_t i=0;i<=N;++i){
double x=beam.node_x(i);
double w=std::cos(w0*x);
double w1=-w0*std::sin(w0*x);
double w2=-w0*w0*std::cos(w0*x);
double w3= w0*w0*w0*std::sin(w0*x);
double w4= w0*w0*w0*w0*std::cos(w0*x);
double EIx=E0*(1+b*x), EIp=E0*b;
double M2=2*EIp*w3+EIx*w4; // EI''=0
ei[i]=EIx; kk[i]=k0; qq[i]=M2+k0*w; ex[i]=w;
(void)w1;
}
beam.set_ei_nodal(ei); beam.set_foundation_nodal(kk); beam.set_q_nodal(qq);
// free ends: set applied M and V consistent with manufactured solution
double x0=0.0,xL=L;
auto setend=[&](bool left,double x){
double EIx=E0*(1+b*x),EIp=E0*b;
double w2=-w0*w0*std::cos(w0*x);
double w3= w0*w0*w0*std::sin(w0*x);
double M=EIx*w2;
double V=EIp*w2+EIx*w3; // (EI w'')'
beam.set_end_moment(M,left);
// applied_end_shear(left) must equal V -> -P=V -> P=-V
beam.set_end_load(-V,left);
};
setend(true,x0); setend(false,xL);
BeamSolver s(beam); auto r=s.solve();
err=maxerr(r.deflection,ex);
}
int main(){
caseA();
caseB();
printf("CaseC clamped-clamped varEI+foundation convergence:\n");
double e1,e2,e3;
caseC(40,e1); caseC(80,e2); caseC(160,e3);
printf(" N=40 err=%.3e N=80 err=%.3e (rate %.2f) N=160 err=%.3e (rate %.2f)\n",
e1,e2,std::log2(e1/e2),e3,std::log2(e2/e3));
printf("CaseD free-free varEI+foundation convergence:\n");
double d1,d2,d3;
caseD(40,d1); caseD(80,d2); caseD(160,d3);
printf(" N=40 err=%.3e N=80 err=%.3e (rate %.2f) N=160 err=%.3e (rate %.2f)\n",
d1,d2,std::log2(d1/d2),d3,std::log2(d2/d3));
return 0;
}
cd /opt/beamfd && g++ -std=c++17 -O2 -Iinclude /tmp/verify.cpp src/beam.cpp src/solver.cpp src/linalg.cpp -o /tmp/verify && /tmp/verify
CaseA cantilever cubic: maxerr=1.000e-02 (w(L)=15.990000 exp 16.000000) CaseB pinned UDL: wmax=2.106475e-02 exact=2.109375e-02 rel=1.375e-03 CaseC clamped-clamped varEI+foundation convergence: N=40 err=9.522e-03 N=80 err=2.427e-03 (rate 1.97) N=160 err=6.124e-04 (rate 1.99) CaseD free-free varEI+foundation convergence: N=40 err=1.207e+00 N=80 err=5.486e-01 (rate 1.14) N=160 err=2.619e-01 (rate 1.07) [stdout] CaseA cantilever cubic: maxerr=1.000e-02 (w(L)=15.990000 exp 16.000000) CaseB pinned UDL: wmax=2.106475e-02 exact=2.109375e-02 rel=1.375e-03 CaseC clamped-clamped varEI+foundation convergence: N=40 err=9.522e-03 N=80 err=2.427e-03 (rate 1.97) N=160 err=6.124e-04 (rate 1.99) CaseD free-free varEI+foundation convergence: N=40 err=1.207e+00 N=80 err=5.486e-01 (rate 1.14) N=160 err=2.619e-01 (rate 1.07)
/opt/beamfd/src/solver.cpp
99 // ---- Interior equilibrium: (EI w'')'' + k w = q (rows i = 2 .. N-2) ------
100 for (std::size_t i = 2; i + 2 < n; ++i) {
101 const double eim = beam_.EI_at(i - 1);
102 const double ei = beam_.EI_at(i);
103 const double eip = beam_.EI_at(i + 1);
104 A(i, i - 2) += eim / h4;
105 A(i, i - 1) += (-2.0 * eim - 2.0 * ei) / h4;
106 A(i, i) += (eim + 4.0 * ei + eip) / h4 + beam_.k_at(i);
107 A(i, i + 1) += (-2.0 * ei - 2.0 * eip) / h4;
108 A(i, i + 2) += eip / h4;
109 sys.b[i] = beam_.q_at(i);
110 }
111
112 // ---- Boundary rows -------------------------------------------------------
113 // Each end contributes its two conditions. The left end fills rows {0, 1}
114 // (which reference low-index nodes), the right end fills rows {n-1, n-2}.
115
116 // w = 0 at node e.
117 auto set_deflection_zero = [&](std::size_t r, std::size_t e) {
118 A(r, e) += 1.0;
119 sys.b[r] = 0.0;
120 };
121 // w' = 0 at an end, via a second-order one-sided first difference.
122 auto set_slope_zero = [&](std::size_t r, bool left) {
123 if (left) {
124 A(r, 0) += -3.0;
125 A(r, 1) += 4.0;
126 A(r, 2) += -1.0;
127 } else {
128 A(r, n - 1) += 3.0;
129 A(r, n - 2) += -4.0;
130 A(r, n - 3) += 1.0;
131 }
132 sys.b[r] = 0.0;
133 };
134 // EI w'' = M at node e (M = 0 for a pinned support).
135 auto set_moment = [&](std::size_t r, std::size_t e, double M) {
136 add_curvature(r, e, 1.0);
137 sys.b[r] = M;
138 };
139 // (EI w'')' = V at an end, via a one-sided first difference of the nodal
140 // moment field, matching shear_field.
141 auto set_shear = [&](std::size_t r, bool left, double V) {
142 const double inv2h = 1.0 / (2.0 * h);
143 if (left) {
144 add_curvature(r, 0, -3.0 * inv2h);
145 add_curvature(r, 1, 4.0 * inv2h);
146 add_curvature(r, 2, -1.0 * inv2h);
147 } else {
148 add_curvature(r, n - 1, 3.0 * inv2h);
149 add_curvature(r, n - 2, -4.0 * inv2h);
150 add_curvature(r, n - 3, 1.0 * inv2h);
151 }
152 sys.b[r] = V;
153 };
154
155 auto apply_end = [&](bool left) {
156 const Support s = left ? beam_.left() : beam_.right();
157 const std::size_t e = left ? 0u : n - 1;
158 const std::size_t r_first = left ? 0u : n - 1; // "primary" condition
159 const std::size_t r_second = left ? 1u : n - 2; // "secondary" condition
160 switch (s) {
161 case Support::Clamped:
162 set_deflection_zero(r_first, e);
163 set_slope_zero(r_second, left);
164 break;
165 case Support::Pinned:
166 set_deflection_zero(r_first, e);
167 set_moment(r_second, e, 0.0);
168 break;
169 case Support::Free:
170 set_moment(r_first, e, beam_.end_moment(left));
171 set_shear(r_second, left, beam_.applied_end_shear(left));
172 break;
173 }
174 };
175 apply_end(true);
176 apply_end(false);
177
178 return sys;
179 }
180
181 BeamResult BeamSolver::solve() const {
182 beam_.validate(); // propagate mechanism / ill-posed configuration errors
183 const BeamSystem sys = assemble_system();
184 BeamResult result;
185 // Propagates std::runtime_error from solve_general if A is singular.
186 result.deflection = solve_general(sys.A, sys.b);
187 return result;
188 }
189
190 // ----------------------------------------------------------------------------
191 // Recovery helpers (already implemented; do NOT modify). These differentiate a
192 // given deflection field so callers can recompute physical quantities (bending
193 // moment, shear) directly from a candidate result, honoring a variable EI(x).
/opt/beamfd/src/solver.cpp
// ---- Boundary rows -------------------------------------------------------
// Each end contributes its two conditions. The left end fills rows {0, 1}
// (which reference low-index nodes), the right end fills rows {n-1, n-2}.
// w = 0 at node e.
auto set_deflection_zero = [&](std::size_t r, std::size_t e) {
A(r, e) += 1.0;
sys.b[r] = 0.0;
};
// w' = 0 at an end, via a second-order one-sided first difference.
auto set_slope_zero = [&](std::size_t r, bool left) {
if (left) {
A(r, 0) += -3.0;
A(r, 1) += 4.0;
A(r, 2) += -1.0;
} else {
A(r, n - 1) += 3.0;
A(r, n - 2) += -4.0;
A(r, n - 3) += 1.0;
}
sys.b[r] = 0.0;
};
// EI w'' = M at node e (M = 0 for a pinned support).
auto set_moment = [&](std::size_t r, std::size_t e, double M) {
add_curvature(r, e, 1.0);
sys.b[r] = M;
};
// (EI w'')' = V at an end, via a one-sided first difference of the nodal
// moment field, matching shear_field.
auto set_shear = [&](std::size_t r, bool left, double V) {
const double inv2h = 1.0 / (2.0 * h);
if (left) {
add_curvature(r, 0, -3.0 * inv2h);
add_curvature(r, 1, 4.0 * inv2h);
add_curvature(r, 2, -1.0 * inv2h);
} else {
add_curvature(r, n - 1, 3.0 * inv2h);
add_curvature(r, n - 2, -4.0 * inv2h);
add_curvature(r, n - 3, 1.0 * inv2h);
}
sys.b[r] = V;
};
auto apply_end = [&](bool left) {
const Support s = left ? beam_.left() : beam_.right();
const std::size_t e = left ? 0u : n - 1;
const std::size_t r_first = left ? 0u : n - 1; // "primary" condition
const std::size_t r_second = left ? 1u : n - 2; // "secondary" condition
switch (s) {
case Support::Clamped:
set_deflection_zero(r_first, e);
set_slope_zero(r_second, left);
break;
case Support::Pinned:
set_deflection_zero(r_first, e);
set_moment(r_second, e, 0.0);
break;
case Support::Free:
set_moment(r_first, e, beam_.end_moment(left));
set_shear(r_second, left, beam_.applied_end_shear(left));
break;
}
};
apply_end(true);
apply_end(false);
return sys; // ---- Boundary rows -------------------------------------------------------
// Each end supplies its two conditions. The left end fills rows {0, 1}, the
// right end fills rows {n-1, n-2}.
//
// Essential conditions (clamped/pinned) are imposed directly. Natural
// conditions (free end: prescribed moment M and shear V) are imposed by the
// ghost-node method: the conservative PDE stencil is applied AT the two
// boundary nodes, and the two ghost values are eliminated using the natural
// conditions written as *central* differences about the end node
//
// M_end = EI_end (w_{g} - 2 w_end + w_in) / h^2 (moment, central)
// V_end = (M_{in} - M_{g}) / (2h) (shear, central)
//
// Substituting collapses the ghost (and ghost-EI) terms exactly, leaving a
// second-order accurate equation in the physical nodes only. This keeps the
// whole scheme O(h^2) at a free end, where a one-sided nesting of the shear
// operator would degrade to O(h).
// w = 0 at node e.
auto set_deflection_zero = [&](std::size_t r, std::size_t e) {
A(r, e) += 1.0;
sys.b[r] = 0.0;
};
// w' = 0 at an end, via a second-order one-sided first difference.
auto set_slope_zero = [&](std::size_t r, bool left) {
if (left) {
A(r, 0) += -3.0;
A(r, 1) += 4.0;
A(r, 2) += -1.0;
} else {
A(r, n - 1) += 3.0;
A(r, n - 2) += -4.0;
A(r, n - 3) += 1.0;
}
sys.b[r] = 0.0;
};
// EI w'' = 0 at node e (pinned moment release), one-sided second difference.
auto set_moment_release = [&](std::size_t r, std::size_t e) {
add_curvature(r, e, 1.0);
sys.b[r] = 0.0;
};
// Free LEFT end: ghost-eliminated PDE at nodes 0 and 1.
auto set_free_left = [&]() {
const double M = beam_.end_moment(true);
const double V = beam_.applied_end_shear(true);
const double ei1 = beam_.EI_at(1);
const double ei2 = beam_.EI_at(2);
// Node 0: 2 EI_1 (w0 - 2 w1 + w2)/h^4 + k0 w0 = q0 + 2 M/h^2 + 2 V/h
A(0, 0) += 2.0 * ei1 / h4 + beam_.k_at(0);
A(0, 1) += -4.0 * ei1 / h4;
A(0, 2) += 2.0 * ei1 / h4;
sys.b[0] = beam_.q_at(0) + 2.0 * M / h2 + 2.0 * V / h;
// Node 1: [-2 EI_1(w0-2w1+w2) + EI_2(w1-2w2+w3)]/h^4 + k1 w1 = q1 - M/h^2
A(1, 0) += -2.0 * ei1 / h4;
A(1, 1) += (4.0 * ei1 + ei2) / h4 + beam_.k_at(1);
A(1, 2) += (-2.0 * ei1 - 2.0 * ei2) / h4;
A(1, 3) += ei2 / h4;
sys.b[1] = beam_.q_at(1) - M / h2;
};
// Free RIGHT end: mirror image at nodes N and N-1.
auto set_free_right = [&]() {
const std::size_t N = n - 1;
const double M = beam_.end_moment(false);
const double V = beam_.applied_end_shear(false);
const double eiN1 = beam_.EI_at(N - 1);
const double eiN2 = beam_.EI_at(N - 2);
// Node N: 2 EI_{N-1}(w_{N-2}-2w_{N-1}+w_N)/h^4 + kN wN = qN + 2 M/h^2 - 2 V/h
A(N, N - 2) += 2.0 * eiN1 / h4;
A(N, N - 1) += -4.0 * eiN1 / h4;
A(N, N) += 2.0 * eiN1 / h4 + beam_.k_at(N);
sys.b[N] = beam_.q_at(N) + 2.0 * M / h2 - 2.0 * V / h;
// Node N-1: [EI_{N-2}(w_{N-3}-2w_{N-2}+w_{N-1})
// - 2 EI_{N-1}(w_{N-2}-2w_{N-1}+w_N)]/h^4 + k_{N-1} w_{N-1}
// = q_{N-1} - M/h^2
A(N - 1, N - 3) += eiN2 / h4;
A(N - 1, N - 2) += (-2.0 * eiN2 - 2.0 * eiN1) / h4;
A(N - 1, N - 1) += (eiN2 + 4.0 * eiN1) / h4 + beam_.k_at(N - 1);
A(N - 1, N) += -2.0 * eiN1 / h4;
sys.b[N - 1] = beam_.q_at(N - 1) - M / h2;
};
auto apply_end = [&](bool left) {
const Support s = left ? beam_.left() : beam_.right();
const std::size_t e = left ? 0u : n - 1;
const std::size_t r_first = left ? 0u : n - 1; // w = 0 row
const std::size_t r_second = left ? 1u : n - 2; // slope / moment row
switch (s) {
case Support::Clamped:
set_deflection_zero(r_first, e);
set_slope_zero(r_second, left);
break;
case Support::Pinned:
set_deflection_zero(r_first, e);
set_moment_release(r_second, e);
break;
case Support::Free:
if (left) {
set_free_left();
} else {
set_free_right();
}
break;
}
};
apply_end(true);
apply_end(false);
return sys;g++ -std=c++17 -O2 -Iinclude /tmp/verify.cpp src/beam.cpp src/solver.cpp src/linalg.cpp -o /tmp/verify && /tmp/verify
CaseA cantilever cubic: maxerr=1.000e-02 (w(L)=15.990000 exp 16.000000) CaseB pinned UDL: wmax=2.106475e-02 exact=2.109375e-02 rel=1.375e-03 CaseC clamped-clamped varEI+foundation convergence: N=40 err=9.522e-03 N=80 err=2.427e-03 (rate 1.97) N=160 err=6.124e-04 (rate 1.99) CaseD free-free varEI+foundation convergence: N=40 err=3.005e-03 N=80 err=7.517e-04 (rate 2.00) N=160 err=1.880e-04 (rate 2.00) [stdout] CaseA cantilever cubic: maxerr=1.000e-02 (w(L)=15.990000 exp 16.000000) CaseB pinned UDL: wmax=2.106475e-02 exact=2.109375e-02 rel=1.375e-03 CaseC clamped-clamped varEI+foundation convergence: N=40 err=9.522e-03 N=80 err=2.427e-03 (rate 1.97) N=160 err=6.124e-04 (rate 1.99) CaseD free-free varEI+foundation convergence: N=40 err=3.005e-03 N=80 err=7.517e-04 (rate 2.00) N=160 err=1.880e-04 (rate 2.00)
#include <cmath>
#include <cstdio>
#include <functional>
#include <vector>
#include "beamfd/beam.hpp"
#include "beamfd/solver.hpp"
using namespace beamfd;
struct WD { double w,w1,w2,w3,w4; };
static double maxerr(const std::vector<double>&a,const std::vector<double>&b){
double m=0; for(size_t i=0;i<a.size();++i) m=std::max(m,std::fabs(a[i]-b[i])); return m;}
// run one manufactured case; returns max error
double run(const char* name, Support L, Support R, std::function<WD(double)> fw,
double Lspan, double E0, double bslope, double k0, size_t N, bool print=false){
Beam beam(Lspan,E0,N,L,R);
std::vector<double> ei(N+1),kk(N+1),qq(N+1),ex(N+1);
for(size_t i=0;i<=N;++i){
double x=beam.node_x(i); WD d=fw(x);
double EIx=E0*(1+bslope*x), EIp=E0*bslope; // EI''=0
double q = EIp*d.w3*2.0 + EIx*d.w4 + k0*d.w; // 2EI'w'''+EI w''''+k w
ei[i]=EIx; kk[i]=k0; qq[i]=q; ex[i]=d.w;
}
beam.set_ei_nodal(ei);
if(k0>0) beam.set_foundation_nodal(kk);
beam.set_q_nodal(qq);
auto endset=[&](bool left,double x){
WD d=fw(x); double EIx=E0*(1+bslope*x),EIp=E0*bslope;
double M=EIx*d.w2;
double V=EIp*d.w2+EIx*d.w3; // (EI w'')'
beam.set_end_moment(M,left);
beam.set_end_load(-V,left); // applied_end_shear=-P must equal V
};
if(L==Support::Free) endset(true,0.0);
if(R==Support::Free) endset(false,Lspan);
BeamSolver s(beam); auto r=s.solve();
double e=maxerr(r.deflection,ex);
if(print) printf(" %s err=%.4e\n",name,e);
return e;
}
void conv(const char* name, Support L, Support R, std::function<WD(double)> fw,
double Lspan,double E0,double b,double k0){
double e1=run(name,L,R,fw,Lspan,E0,b,k0,40);
double e2=run(name,L,R,fw,Lspan,E0,b,k0,80);
double e3=run(name,L,R,fw,Lspan,E0,b,k0,160);
printf("%-16s N40=%.3e N80=%.3e(%.2f) N160=%.3e(%.2f)\n",
name,e1,e2,std::log2(e1/e2),e3,std::log2(e2/e3));
}
int main(){
double Ls=2.0,E0=700.0,b=0.4/2.0,k0=250.0;
double pi=std::acos(-1.0);
// clamped-clamped: w=x^2(L-x)^2
conv("clamp-clamp",Support::Clamped,Support::Clamped,[&](double x){
double L=Ls; return WD{ x*x*(L-x)*(L-x),
2*x*(L-x)*(L-x)-2*x*x*(L-x),
2*L*L-12*L*x+12*x*x, -12*L+24*x, 24.0};},Ls,E0,b,k0);
// pinned-pinned: w=sin(pi x/L)
conv("pin-pin",Support::Pinned,Support::Pinned,[&](double x){
double a=pi/Ls; return WD{std::sin(a*x),a*std::cos(a*x),-a*a*std::sin(a*x),
-a*a*a*std::cos(a*x),a*a*a*a*std::sin(a*x)};},Ls,E0,b,k0);
// free-free: w=cos(pi x/L)
conv("free-free",Support::Free,Support::Free,[&](double x){
double a=pi/Ls; return WD{std::cos(a*x),-a*std::sin(a*x),-a*a*std::cos(a*x),
a*a*a*std::sin(a*x),a*a*a*a*std::cos(a*x)};},Ls,E0,b,k0);
// clamped-free (cantilever): w=1-cos(pi x/(2L)) -> w(0)=0,w'(0)=0
conv("clamp-free",Support::Clamped,Support::Free,[&](double x){
double a=pi/(2*Ls); return WD{1-std::cos(a*x),a*std::sin(a*x),a*a*std::cos(a*x),
-a*a*a*std::sin(a*x),-a*a*a*a*std::cos(a*x)};},Ls,E0,b,k0);
// pinned-free: w=sin(pi x/(2L)) -> w(0)=0,w''(0)=0
conv("pin-free",Support::Pinned,Support::Free,[&](double x){
double a=pi/(2*Ls); return WD{std::sin(a*x),a*std::cos(a*x),-a*a*std::sin(a*x),
-a*a*a*std::cos(a*x),a*a*a*a*std::sin(a*x)};},Ls,E0,b,k0);
// clamped-pinned: quartic w with w(0)=0,w'(0)=0,w(L)=0,w''(L)=0
// pick w=x^2(L-x)(c) ... use w = x^2 (L-x)(a x + c). Need w''(L)=0.
// Simpler: w = x^2 (L - x)^2 + d * x^2 (L-x). Let's just use cubic-ish:
// w = x^2(L-x): w(0)=0,w'(0)=0,w(L)=0, w''=2L-6x, w''(L)=2L-6L=-4L != 0. no.
// Use w = x^2 (L-x)^2 (satisfies clamped both & also w''(L)=2L^2-12L^2+12L^2=2L^2!=0)
// -> not pinned at right. Build explicit quartic satisfying 4 conds:
// w=a x^2 + ... too messy; use w = sin via shift: w=x^2*(L-x)^2 won't pin.
// Construct: w(x)=x^2 (L-x) (L - x/... ) -- do numerically below.
conv("clamp-pin",Support::Clamped,Support::Pinned,[&](double x){
// w = x^2(L-x)^2*(... )? Need w''(L)=0. Use w = x^2(L-x)^2 - (1/6)*?
// Instead use w = x^2*(L-x)* (3L-2x)/(something)... derive quartic:
// Let w=A x^2 + B x^3 + C x^4 with w(L)=0,w''(L)=0 (w(0)=0,w'(0)=0 auto).
// w(L)=A L^2+B L^3+C L^4=0 ; w''=2A+6Bx+12Cx^2; w''(L)=2A+6BL+12CL^2=0.
double L=Ls;
double C=1.0; // free scale
// solve A,B from: A L^2 + B L^3 = -C L^4 ; 2A + 6 B L = -12 C L^2
// => A + B L = -C L^2 ; A + 3 B L = -6 C L^2 ; subtract: 2BL=-5CL^2 => B=-2.5 C L
double B=-2.5*C*L;
double A=-C*L*L - B*L; // = -C L^2 +2.5 C L^2 = 1.5 C L^2
double w=A*x*x+B*x*x*x+C*x*x*x*x;
double w1=2*A*x+3*B*x*x+4*C*x*x*x;
double w2=2*A+6*B*x+12*C*x*x;
double w3=6*B+24*C*x;
double w4=24*C;
return WD{w,w1,w2,w3,w4};},Ls,E0,b,k0);
return 0;
}
g++ -std=c++17 -O2 -Iinclude /tmp/verify2.cpp src/beam.cpp src/solver.cpp src/linalg.cpp -o /tmp/verify2 && /tmp/verify2
clamp-clamp N40=9.523e-03 N80=2.427e-03(1.97) N160=6.125e-04(1.99) pin-pin N40=4.112e-04 N80=1.737e-04(1.24) N160=5.230e-05(1.73) free-free N40=1.499e-02 N80=3.750e-03(2.00) N160=9.377e-04(2.00) clamp-free N40=2.903e-04 N80=7.453e-05(1.96) N160=1.888e-05(1.98) pin-free N40=4.311e-04 N80=1.169e-04(1.88) N160=3.040e-05(1.94) clamp-pin N40=1.650e-02 N80=4.158e-03(1.99) N160=1.044e-03(1.99) [stdout] clamp-clamp N40=9.523e-03 N80=2.427e-03(1.97) N160=6.125e-04(1.99) pin-pin N40=4.112e-04 N80=1.737e-04(1.24) N160=5.230e-05(1.73) free-free N40=1.499e-02 N80=3.750e-03(2.00) N160=9.377e-04(2.00) clamp-free N40=2.903e-04 N80=7.453e-05(1.96) N160=1.888e-05(1.98) pin-free N40=4.311e-04 N80=1.169e-04(1.88) N160=3.040e-05(1.94) clamp-pin N40=1.650e-02 N80=4.158e-03(1.99) N160=1.044e-03(1.99)
cat > /tmp/pp.cpp <<'EOF'
#include <cmath>
#include <cstdio>
#include <vector>
#include "beamfd/beam.hpp"
#include "beamfd/solver.hpp"
using namespace beamfd;
double run(size_t N){
double Ls=2.0,E0=700.0,b=0.2,k0=250.0,pi=std::acos(-1.0),a=pi/Ls;
Beam beam(Ls,E0,N,Support::Pinned,Support::Pinned);
std::vector<double> ei(N+1),kk(N+1),qq(N+1),ex(N+1);
for(size_t i=0;i<=N;++i){double x=beam.node_x(i);
double w=std::sin(a*x),w3=-a*a*a*std::cos(a*x),w4=a*a*a*a*std::sin(a*x);
double EIx=E0*(1+b*x),EIp=E0*b;
ei[i]=EIx;kk[i]=k0;qq[i]=2*EIp*w3+EIx*w4+k0*w;ex[i]=w;}
beam.set_ei_nodal(ei);beam.set_foundation_nodal(kk);beam.set_q_nodal(qq);
BeamSolver s(beam);auto r=s.solve();
double m=0;for(size_t i=0;i<=N;++i)m=std::max(m,std::fabs(r.deflection[i]-ex[i]));return m;}
int main(){double p=0;for(size_t N=40;N<=2560;N*=2){double e=run(N);
printf("N=%4zu err=%.5e rate=%.3f\n",N,e,p>0?std::log2(p/e):0);p=e;}}
EOF
g++ -std=c++17 -O2 -Iinclude /tmp/pp.cpp src/beam.cpp src/solver.cpp src/linalg.cpp -o /tmp/pp && /tmp/ppN= 40 err=4.11225e-04 rate=0.000 N= 80 err=1.73664e-04 rate=1.244 N= 160 err=5.22978e-05 rate=1.731 N= 320 err=1.41819e-05 rate=1.883 N= 640 err=3.65722e-06 rate=1.955 N=1280 err=4.10487e-07 rate=3.155 N=2560 err=1.90304e-05 rate=-5.535 [stdout] N= 40 err=4.11225e-04 rate=0.000 N= 80 err=1.73664e-04 rate=1.244 N= 160 err=5.22978e-05 rate=1.731 N= 320 err=1.41819e-05 rate=1.883 N= 640 err=3.65722e-06 rate=1.955 N=1280 err=4.10487e-07 rate=3.155 N=2560 err=1.90304e-05 rate=-5.535
/opt/beamfd/src/solver.cpp
112 // ---- Boundary rows -------------------------------------------------------
113 // Each end supplies its two conditions. The left end fills rows {0, 1}, the
114 // right end fills rows {n-1, n-2}.
115 //
116 // Essential conditions (clamped/pinned) are imposed directly. Natural
117 // conditions (free end: prescribed moment M and shear V) are imposed by the
118 // ghost-node method: the conservative PDE stencil is applied AT the two
119 // boundary nodes, and the two ghost values are eliminated using the natural
120 // conditions written as *central* differences about the end node
121 //
122 // M_end = EI_end (w_{g} - 2 w_end + w_in) / h^2 (moment, central)
123 // V_end = (M_{in} - M_{g}) / (2h) (shear, central)
124 //
125 // Substituting collapses the ghost (and ghost-EI) terms exactly, leaving a
126 // second-order accurate equation in the physical nodes only. This keeps the
127 // whole scheme O(h^2) at a free end, where a one-sided nesting of the shear
128 // operator would degrade to O(h).
129
130 // w = 0 at node e.
131 auto set_deflection_zero = [&](std::size_t r, std::size_t e) {
132 A(r, e) += 1.0;
133 sys.b[r] = 0.0;
134 };
135 // w' = 0 at an end, via a second-order one-sided first difference.
136 auto set_slope_zero = [&](std::size_t r, bool left) {
137 if (left) {
138 A(r, 0) += -3.0;
139 A(r, 1) += 4.0;
140 A(r, 2) += -1.0;
141 } else {
142 A(r, n - 1) += 3.0;
143 A(r, n - 2) += -4.0;
144 A(r, n - 3) += 1.0;
145 }
146 sys.b[r] = 0.0;
147 };
148 // EI w'' = 0 at node e (pinned moment release), one-sided second difference.
149 auto set_moment_release = [&](std::size_t r, std::size_t e) {
150 add_curvature(r, e, 1.0);
151 sys.b[r] = 0.0;
152 };
153
154 // Free LEFT end: ghost-eliminated PDE at nodes 0 and 1.
155 auto set_free_left = [&]() {
156 const double M = beam_.end_moment(true);
157 const double V = beam_.applied_end_shear(true);
158 const double ei1 = beam_.EI_at(1);
159 const double ei2 = beam_.EI_at(2);
160 // Node 0: 2 EI_1 (w0 - 2 w1 + w2)/h^4 + k0 w0 = q0 + 2 M/h^2 + 2 V/h
161 A(0, 0) += 2.0 * ei1 / h4 + beam_.k_at(0);
162 A(0, 1) += -4.0 * ei1 / h4;
163 A(0, 2) += 2.0 * ei1 / h4;
164 sys.b[0] = beam_.q_at(0) + 2.0 * M / h2 + 2.0 * V / h;
165 // Node 1: [-2 EI_1(w0-2w1+w2) + EI_2(w1-2w2+w3)]/h^4 + k1 w1 = q1 - M/h^2
166 A(1, 0) += -2.0 * ei1 / h4;
167 A(1, 1) += (4.0 * ei1 + ei2) / h4 + beam_.k_at(1);
168 A(1, 2) += (-2.0 * ei1 - 2.0 * ei2) / h4;
169 A(1, 3) += ei2 / h4;
170 sys.b[1] = beam_.q_at(1) - M / h2;
171 };
172 // Free RIGHT end: mirror image at nodes N and N-1.
173 auto set_free_right = [&]() {
174 const std::size_t N = n - 1;
175 const double M = beam_.end_moment(false);
176 const double V = beam_.applied_end_shear(false);
177 const double eiN1 = beam_.EI_at(N - 1);
178 const double eiN2 = beam_.EI_at(N - 2);
179 // Node N: 2 EI_{N-1}(w_{N-2}-2w_{N-1}+w_N)/h^4 + kN wN = qN + 2 M/h^2 - 2 V/h
180 A(N, N - 2) += 2.0 * eiN1 / h4;
181 A(N, N - 1) += -4.0 * eiN1 / h4;
182 A(N, N) += 2.0 * eiN1 / h4 + beam_.k_at(N);
183 sys.b[N] = beam_.q_at(N) + 2.0 * M / h2 - 2.0 * V / h;
184 // Node N-1: [EI_{N-2}(w_{N-3}-2w_{N-2}+w_{N-1})
185 // - 2 EI_{N-1}(w_{N-2}-2w_{N-1}+w_N)]/h^4 + k_{N-1} w_{N-1}
186 // = q_{N-1} - M/h^2
187 A(N - 1, N - 3) += eiN2 / h4;
188 A(N - 1, N - 2) += (-2.0 * eiN2 - 2.0 * eiN1) / h4;
189 A(N - 1, N - 1) += (eiN2 + 4.0 * eiN1) / h4 + beam_.k_at(N - 1);
190 A(N - 1, N) += -2.0 * eiN1 / h4;
191 sys.b[N - 1] = beam_.q_at(N - 1) - M / h2;
192 };
193
194 auto apply_end = [&](bool left) {
195 const Support s = left ? beam_.left() : beam_.right();
196 const std::size_t e = left ? 0u : n - 1;
197 const std::size_t r_first = left ? 0u : n - 1; // w = 0 row
198 const std::size_t r_second = left ? 1u : n - 2; // slope / moment row
199 switch (s) {
200 case Support::Clamped:
201 set_deflection_zero(r_first, e);
202 set_slope_zero(r_second, left);
203 break;
204 case Support::Pinned:
205 set_deflection_zero(r_first, e);
206 set_moment_release(r_second, e);
207 break;
208 case Support::Free:
209 if (left) {
210 set_free_left();
211 } else {
212 set_free_right();
213 }
214 break;
215 }
216 };
217 apply_end(true);
218 apply_end(false);
219
220 return sys;
221 }
222
223 BeamResult BeamSolver::solve() const {
224 beam_.validate(); // propagate mechanism / ill-posed configuration errors
225 const BeamSystem sys = assemble_system();
226 BeamResult result;
227 // Propagates std::runtime_error from solve_general if A is singular.
228 result.deflection = solve_general(sys.A, sys.b);
229 return result;
230 }
231
/opt/beamfd/src/solver.cpp
// EI w'' = 0 at node e (pinned moment release), one-sided second difference.
auto set_moment_release = [&](std::size_t r, std::size_t e) {
add_curvature(r, e, 1.0);
sys.b[r] = 0.0;
};
// Free LEFT end: ghost-eliminated PDE at nodes 0 and 1.
auto set_free_left = [&]() {
const double M = beam_.end_moment(true);
const double V = beam_.applied_end_shear(true);
const double ei1 = beam_.EI_at(1);
const double ei2 = beam_.EI_at(2);
// Node 0: 2 EI_1 (w0 - 2 w1 + w2)/h^4 + k0 w0 = q0 + 2 M/h^2 + 2 V/h
A(0, 0) += 2.0 * ei1 / h4 + beam_.k_at(0);
A(0, 1) += -4.0 * ei1 / h4;
A(0, 2) += 2.0 * ei1 / h4;
sys.b[0] = beam_.q_at(0) + 2.0 * M / h2 + 2.0 * V / h;
// Node 1: [-2 EI_1(w0-2w1+w2) + EI_2(w1-2w2+w3)]/h^4 + k1 w1 = q1 - M/h^2
A(1, 0) += -2.0 * ei1 / h4;
A(1, 1) += (4.0 * ei1 + ei2) / h4 + beam_.k_at(1);
A(1, 2) += (-2.0 * ei1 - 2.0 * ei2) / h4;
A(1, 3) += ei2 / h4;
sys.b[1] = beam_.q_at(1) - M / h2;
};
// Free RIGHT end: mirror image at nodes N and N-1.
auto set_free_right = [&]() {
const std::size_t N = n - 1;
const double M = beam_.end_moment(false);
const double V = beam_.applied_end_shear(false);
const double eiN1 = beam_.EI_at(N - 1);
const double eiN2 = beam_.EI_at(N - 2);
// Node N: 2 EI_{N-1}(w_{N-2}-2w_{N-1}+w_N)/h^4 + kN wN = qN + 2 M/h^2 - 2 V/h
A(N, N - 2) += 2.0 * eiN1 / h4;
A(N, N - 1) += -4.0 * eiN1 / h4;
A(N, N) += 2.0 * eiN1 / h4 + beam_.k_at(N);
sys.b[N] = beam_.q_at(N) + 2.0 * M / h2 - 2.0 * V / h;
// Node N-1: [EI_{N-2}(w_{N-3}-2w_{N-2}+w_{N-1})
// - 2 EI_{N-1}(w_{N-2}-2w_{N-1}+w_N)]/h^4 + k_{N-1} w_{N-1}
// = q_{N-1} - M/h^2
A(N - 1, N - 3) += eiN2 / h4;
A(N - 1, N - 2) += (-2.0 * eiN2 - 2.0 * eiN1) / h4;
A(N - 1, N - 1) += (eiN2 + 4.0 * eiN1) / h4 + beam_.k_at(N - 1);
A(N - 1, N) += -2.0 * eiN1 / h4;
sys.b[N - 1] = beam_.q_at(N - 1) - M / h2;
};
auto apply_end = [&](bool left) {
const Support s = left ? beam_.left() : beam_.right();
const std::size_t e = left ? 0u : n - 1;
const std::size_t r_first = left ? 0u : n - 1; // w = 0 row
const std::size_t r_second = left ? 1u : n - 2; // slope / moment row
switch (s) {
case Support::Clamped:
set_deflection_zero(r_first, e);
set_slope_zero(r_second, left);
break;
case Support::Pinned:
set_deflection_zero(r_first, e);
set_moment_release(r_second, e);
break;
case Support::Free:
if (left) {
set_free_left();
} else {
set_free_right();
}
break;
}
}; // Near-boundary PDE row at the first interior node (node 1 on the left,
// node N-1 on the right). It is the conservative stencil with the end-node
// ghost eliminated through the moment condition EI_end w''_end = M (central),
// which folds the end curvature into the right-hand side as -M/h^2. Shared by
// free ends (M = applied end moment) and pinned ends (moment release, M = 0).
auto set_ghost_pde_row = [&](bool left, double M) {
if (left) {
const double ei1 = beam_.EI_at(1);
const double ei2 = beam_.EI_at(2);
A(1, 0) += -2.0 * ei1 / h4;
A(1, 1) += (4.0 * ei1 + ei2) / h4 + beam_.k_at(1);
A(1, 2) += (-2.0 * ei1 - 2.0 * ei2) / h4;
A(1, 3) += ei2 / h4;
sys.b[1] = beam_.q_at(1) - M / h2;
} else {
const std::size_t N = n - 1;
const double eiN1 = beam_.EI_at(N - 1);
const double eiN2 = beam_.EI_at(N - 2);
A(N - 1, N - 3) += eiN2 / h4;
A(N - 1, N - 2) += (-2.0 * eiN2 - 2.0 * eiN1) / h4;
A(N - 1, N - 1) += (eiN2 + 4.0 * eiN1) / h4 + beam_.k_at(N - 1);
A(N - 1, N) += -2.0 * eiN1 / h4;
sys.b[N - 1] = beam_.q_at(N - 1) - M / h2;
}
};
// End-node PDE row for a free end. The conservative stencil at the end node
// with both ghosts eliminated through the moment (M) and shear (V) conditions
// written as central differences; the ghost-EI terms cancel exactly.
auto set_free_end_row = [&](bool left, double M, double V) {
if (left) {
const double ei1 = beam_.EI_at(1);
// 2 EI_1 (w0 - 2 w1 + w2)/h^4 + k0 w0 = q0 + 2 M/h^2 + 2 V/h
A(0, 0) += 2.0 * ei1 / h4 + beam_.k_at(0);
A(0, 1) += -4.0 * ei1 / h4;
A(0, 2) += 2.0 * ei1 / h4;
sys.b[0] = beam_.q_at(0) + 2.0 * M / h2 + 2.0 * V / h;
} else {
const std::size_t N = n - 1;
const double eiN1 = beam_.EI_at(N - 1);
// 2 EI_{N-1}(w_{N-2}-2w_{N-1}+w_N)/h^4 + kN wN = qN + 2 M/h^2 - 2 V/h
A(N, N - 2) += 2.0 * eiN1 / h4;
A(N, N - 1) += -4.0 * eiN1 / h4;
A(N, N) += 2.0 * eiN1 / h4 + beam_.k_at(N);
sys.b[N] = beam_.q_at(N) + 2.0 * M / h2 - 2.0 * V / h;
}
};
auto apply_end = [&](bool left) {
const Support s = left ? beam_.left() : beam_.right();
const std::size_t e = left ? 0u : n - 1;
const std::size_t r_first = left ? 0u : n - 1; // w = 0 / free end-node row
switch (s) {
case Support::Clamped:
set_deflection_zero(r_first, e);
set_slope_zero(left ? 1u : n - 2, left);
break;
case Support::Pinned:
set_deflection_zero(r_first, e);
set_ghost_pde_row(left, 0.0); // moment release: M = 0
break;
case Support::Free:
set_free_end_row(left, beam_.end_moment(left),
beam_.applied_end_shear(left));
set_ghost_pde_row(left, beam_.end_moment(left));
break;
}
};g++ -std=c++17 -O2 -Iinclude /tmp/verify2.cpp src/beam.cpp src/solver.cpp src/linalg.cpp -o /tmp/verify2 && /tmp/verify2 && echo "--- pin-pin fine ---" && g++ -std=c++17 -O2 -Iinclude /tmp/pp.cpp src/beam.cpp src/solver.cpp src/linalg.cpp -o /tmp/pp && /tmp/pp
clamp-clamp N40=9.523e-03 N80=2.427e-03(1.97) N160=6.125e-04(1.99) pin-pin N40=9.794e-04 N80=2.447e-04(2.00) N160=6.118e-05(2.00) free-free N40=1.499e-02 N80=3.750e-03(2.00) N160=9.377e-04(2.00) clamp-free N40=2.903e-04 N80=7.453e-05(1.96) N160=1.888e-05(1.98) pin-free N40=5.042e-04 N80=1.261e-04(2.00) N160=3.154e-05(2.00) clamp-pin N40=8.164e-03 N80=2.074e-03(1.98) N160=5.225e-04(1.99) --- pin-pin fine --- N= 40 err=9.79359e-04 rate=0.000 N= 80 err=2.44745e-04 rate=2.001 N= 160 err=6.11838e-05 rate=2.000 N= 320 err=1.52936e-05 rate=2.000 N= 640 err=3.78886e-06 rate=2.013 N=1280 err=3.39369e-07 rate=3.481 N=2560 err=7.39310e-06 rate=-4.445 [stdout] clamp-clamp N40=9.523e-03 N80=2.427e-03(1.97) N160=6.125e-04(1.99) pin-pin N40=9.794e-04 N80=2.447e-04(2.00) N160=6.118e-05(2.00) free-free N40=1.499e-02 N80=3.750e-03(2.00) N160=9.377e-04(2.00) clamp-free N40=2.903e-04 N80=7.453e-05(1.96) N160=1.888e-05(1.98) pin-free N40=5.042e-04 N80=1.261e-04(2.00) N160=3.154e-05(2.00) clamp-pin N40=8.164e-03 N80=2.074e-03(1.98) N160=5.225e-04(1.99) --- pin-pin fine --- N= 40 err=9.79359e-04 rate=0.000 N= 80 err=2.44745e-04 rate=2.001 N= 160 err=6.11838e-05 rate=2.000 N= 320 err=1.52936e-05 rate=2.000 N= 640 err=3.78886e-06 rate=2.013 N=1280 err=3.39369e-07 rate=3.481 N=2560 err=7.39310e-06 rate=-4.445
cmake --build build --parallel >/tmp/b.log 2>&1 && echo BUILD_OK && cd build && ctest --output-on-failure 2>&1 | tail -6 && echo "--- demo ---" && (./beamfd ../examples/cantilever_tip.beam 2>&1 | head -20 || ls)
BUILD_OK
Start 1: smoke
1/1 Test #1: smoke ............................ Passed 0.00 sec
100% tests passed, 0 tests failed out of 1
Total Test time (real) = 0.00 sec
--- demo ---
/bin/bash: line 1: ./beamfd: No such file or directory
[stdout]
BUILD_OK
Start 1: smoke
1/1 Test #1: smoke ............................ Passed 0.00 sec
100% tests passed, 0 tests failed out of 1
Total Test time (real) = 0.00 sec
--- demo ---
/bin/bash: line 1: ./beamfd: No such file or directoryfind . -maxdepth 2 -type f -executable | grep -v CMake; echo "--- example file ---"; cat ../examples/cantilever_tip.beam
./tests/test_smoke ./beamfd_cli --- example file --- # Cantilever beam with a downward point load at the free tip. # Clamped (built-in) at x = 0, free at x = L. length 3.0 # span [m] ei 2.0e6 # flexural rigidity EI [N*m^2] cells 200 # number of equal finite-difference intervals # Left end clamped, right end free. support clamped free # Downward point load of 1 kN at the free right tip. endload right 1000.0 [stdout] ./tests/test_smoke ./beamfd_cli --- example file --- # Cantilever beam with a downward point load at the free tip. # Clamped (built-in) at x = 0, free at x = L. length 3.0 # span [m] ei 2.0e6 # flexural rigidity EI [N*m^2] cells 200 # number of equal finite-difference intervals # Left end clamped, right end free. support clamped free # Downward point load of 1 kN at the free right tip. endload right 1000.0
./beamfd_cli ../examples/cantilever_tip.beam 2>&1 | head -8; echo "..."; ./beamfd_cli ../examples/cantilever_tip.beam 2>&1 | tail -4
# beam analysis report length=3.000000e+00 EI(max)=2.000000e+06 nodes=201 dx=1.500000e-02 [fields] (i, x, w[m], M=EI*w''[N*m], V=(EI*w'')'[N]) 0 0.000000e+00 -1.596939e-11 3.000000e+03 -1.000000e+03 1 1.500000e-02 1.678905e-07 2.985000e+03 -1.000000e+03 2 3.000000e-02 6.716094e-07 2.970000e+03 -1.000000e+03 3 4.500000e-02 1.509453e-06 2.955000e+03 -1.000000e+03 ... 197 2.955000e+00 4.398647e-03 4.500000e+01 -1.000000e+03 198 2.970000e+00 4.432391e-03 3.000000e+01 -1.000000e+03 199 2.985000e+00 4.466138e-03 1.500000e+01 -1.000000e+03 200 3.000000e+00 4.499887e-03 0.000000e+00 -1.000000e+03 [stdout] # beam analysis report length=3.000000e+00 EI(max)=2.000000e+06 nodes=201 dx=1.500000e-02 [fields] (i, x, w[m], M=EI*w''[N*m], V=(EI*w'')'[N]) 0 0.000000e+00 -1.596939e-11 3.000000e+03 -1.000000e+03 1 1.500000e-02 1.678905e-07 2.985000e+03 -1.000000e+03 2 3.000000e-02 6.716094e-07 2.970000e+03 -1.000000e+03 3 4.500000e-02 1.509453e-06 2.955000e+03 -1.000000e+03 ... 197 2.955000e+00 4.398647e-03 4.500000e+01 -1.000000e+03 198 2.970000e+00 4.432391e-03 3.000000e+01 -1.000000e+03 199 2.985000e+00 4.466138e-03 1.500000e+01 -1.000000e+03 200 3.000000e+00 4.499887e-03 0.000000e+00 -1.000000e+03
cat > /tmp/edge.cpp <<'EOF'
#include <cstdio>
#include <stdexcept>
#include "beamfd/beam.hpp"
#include "beamfd/solver.hpp"
using namespace beamfd;
int main(){
// free-free, no foundation -> mechanism -> validate throws, propagates from solve()
try{ Beam b(2,1000,20,Support::Free,Support::Free); b.set_distributed_load(1);
BeamSolver(b).solve(); printf("FAIL: no throw\n"); }
catch(const std::exception&e){ printf("OK free-free throws: %s\n",e.what()); }
// free-free WITH foundation -> solvable
try{ Beam b(2,1000,20,Support::Free,Support::Free); b.set_foundation_nodal(std::vector<double>(21,500));
b.set_distributed_load(10); auto r=BeamSolver(b).solve();
printf("OK free-free+foundation solves, w_mid=%.4e\n",r.deflection[10]); }
catch(const std::exception&e){ printf("FAIL: %s\n",e.what()); }
// pinned-free with applied end moment at free end
try{ Beam b(3,2e6,50,Support::Pinned,Support::Free);
b.set_end_moment(500.0,false); auto r=BeamSolver(b).solve();
printf("OK pin-free endmoment, w_tip=%.4e\n",r.deflection[50]); }
catch(const std::exception&e){ printf("FAIL: %s\n",e.what()); }
return 0;
}
EOF
cd /opt/beamfd && g++ -std=c++17 -O2 -Iinclude /tmp/edge.cpp src/beam.cpp src/solver.cpp src/linalg.cpp -o /tmp/edge && /tmp/edgeOK free-free throws: Beam::validate: free-free beam is a mechanism OK free-free+foundation solves, w_mid=2.0000e-02 FAIL: Beam::validate: a free end requires a clamped opposite end [stdout] OK free-free throws: Beam::validate: free-free beam is a mechanism OK free-free+foundation solves, w_mid=2.0000e-02 FAIL: Beam::validate: a free end requires a clamped opposite end
/opt/beamfd/src/solver.cpp
39 // ----------------------------------------------------------------------------
40 // Finite-difference discretization of
41 //
42 // (EI(x) w''(x))'' + k(x) w(x) = q(x)
43 //
44 // on the uniform grid x_i = i*dx, i = 0..N, with n = N+1 nodal unknowns.
45 //
46 // Interior nodes (i = 2 .. N-2) use the conservative ("divergence-form") second-
47 // order 5-point stencil obtained by nesting two central second differences:
48 //
49 // M_i = EI_i (w_{i-1} - 2 w_i + w_{i+1}) / h^2 (curvature)
50 // M''_i = (M_{i-1} - 2 M_i + M_{i+1}) / h^2
51 //
52 // which is exact for cubics and O(h^2) in general. The four remaining equations
53 // (two per end) impose the boundary conditions to the same second order, using
54 // exactly the one-sided difference stencils that the recovery helpers above use
55 // for M = EI w'' and V = (EI w'')', so a recovered field reproduces the imposed
56 // moment/shear. Essential conditions (w, w') are imposed directly; natural ones
57 // (M, V at a free end) match the applied end actions.
58 // ----------------------------------------------------------------------------
59
60 BeamSystem BeamSolver::assemble_system() const {
61 const std::size_t n = beam_.num_nodes();
62 if (n < 5) {
63 // The second-order one-sided boundary stencils (and the shear recovery)
64 // need at least four nodes adjacent to each end; the framework's
65 // shear_field carries the same minimum.
66 throw std::runtime_error("BeamSolver::assemble_system: need at least 5 nodes");
67 }
68 const double h = beam_.dx();
69 const double h2 = h * h;
70 const double h4 = h2 * h2;
71
72 BeamSystem sys;
73 sys.A = DenseMatrix(n, n);
74 sys.b.assign(n, 0.0);
75 DenseMatrix& A = sys.A;
76
77 // Add coef * EI(j) * w''(x_j) (using the same second-difference stencil as
78 // second_derivative: one-sided at the ends, central in the interior) to the
79 // matrix row r. The 1/h^2 of the second difference is folded in here.
80 auto add_curvature = [&](std::size_t r, std::size_t j, double coef) {
81 const double f = coef * beam_.EI_at(j) / h2;
82 if (j == 0) {
83 A(r, 0) += 2.0 * f;
84 A(r, 1) += -5.0 * f;
85 A(r, 2) += 4.0 * f;
86 A(r, 3) += -1.0 * f;
87 } else if (j == n - 1) {
88 A(r, n - 1) += 2.0 * f;
89 A(r, n - 2) += -5.0 * f;
90 A(r, n - 3) += 4.0 * f;
91 A(r, n - 4) += -1.0 * f;
92 } else {
93 A(r, j - 1) += 1.0 * f;
94 A(r, j) += -2.0 * f;
95 A(r, j + 1) += 1.0 * f;
96 }
97 };
98
99 // ---- Interior equilibrium: (EI w'')'' + k w = q (rows i = 2 .. N-2) ------
100 for (std::size_t i = 2; i + 2 < n; ++i) {
101 const double eim = beam_.EI_at(i - 1);
102 const double ei = beam_.EI_at(i);
103 const double eip = beam_.EI_at(i + 1);
104 A(i, i - 2) += eim / h4;
105 A(i, i - 1) += (-2.0 * eim - 2.0 * ei) / h4;
106 A(i, i) += (eim + 4.0 * ei + eip) / h4 + beam_.k_at(i);
107 A(i, i + 1) += (-2.0 * ei - 2.0 * eip) / h4;
108 A(i, i + 2) += eip / h4;
109 sys.b[i] = beam_.q_at(i);
110 }
111
112 // ---- Boundary rows -------------------------------------------------------
113 // Each end supplies its two conditions. The left end fills rows {0, 1}, the
/opt/beamfd/src/solver.cpp
// which is exact for cubics and O(h^2) in general. The four remaining equations // (two per end) impose the boundary conditions to the same second order, using // exactly the one-sided difference stencils that the recovery helpers above use // for M = EI w'' and V = (EI w'')', so a recovered field reproduces the imposed // moment/shear. Essential conditions (w, w') are imposed directly; natural ones // (M, V at a free end) match the applied end actions. // ----------------------------------------------------------------------------
// which is exact for cubics and O(h^2) in general. The four remaining equations // (two per end) impose the boundary conditions to the same second order: // * Essential conditions are imposed directly -- w = 0 (clamped/pinned) and, // for a clamped end, w' = 0 via a one-sided second-order first difference. // * Natural conditions (a free end's moment/shear, and a pinned end's moment // release) are imposed by the ghost-node method: the conservative stencil is // applied at the boundary node(s) and the ghosts are eliminated using the // conditions written as central differences, which preserves O(h^2) where a // one-sided nesting of the shear operator would drop to O(h). // ----------------------------------------------------------------------------
cmake --build build --parallel >/tmp/b.log 2>&1 && echo BUILD_OK && cd build && ctest --output-on-failure 2>&1 | tail -4
BUILD_OK 100% tests passed, 0 tests failed out of 1 Total Test time (real) = 0.00 sec [stdout] BUILD_OK 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.4LkWH2z2RJ/beamfd/build_grade
== Building ==
[ 6%] Building CXX object CMakeFiles/beamfd.dir/src/beam.cpp.o
[ 13%] Building CXX object CMakeFiles/beamfd.dir/src/linalg.cpp.o
[ 20%] Building CXX object CMakeFiles/beamfd.dir/src/io.cpp.o
[ 26%] Building CXX object CMakeFiles/beamfd.dir/src/solver.cpp.o
[ 33%] Linking CXX static library libbeamfd.a
[ 33%] Built target beamfd
[ 40%] Building CXX object tests/CMakeFiles/test_mms.dir/test_mms.cpp.o
[ 46%] Building CXX object CMakeFiles/beamfd_cli.dir/src/main.cpp.o
[ 53%] Building CXX object tests/CMakeFiles/test_soak.dir/test_soak.cpp.o
[ 60%] Building CXX object tests/CMakeFiles/test_convergence.dir/test_convergence.cpp.o
[ 66%] Building CXX object tests/CMakeFiles/test_invariants.dir/test_invariants.cpp.o
[ 73%] Linking CXX executable beamfd_cli
[ 73%] Built target beamfd_cli
[ 80%] Linking CXX executable test_invariants
[ 86%] Linking CXX executable test_mms
[ 93%] Linking CXX executable test_convergence
[ 93%] Built target test_invariants
[100%] Linking CXX executable test_soak
[100%] Built target test_mms
[100%] Built target test_convergence
[100%] Built target test_soak
== Running hidden test suite (ctest) ==
Test project /tmp/tmp.4LkWH2z2RJ/beamfd/build_grade
Start 1: test_mms
1/4 Test #1: test_mms ......................... Passed 0.00 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 ..................***Failed 0.00 sec
[ PASS ] invariant_equilibrium_residual
[ FAIL ] invariant_symmetry: symmetry at node 0
[ PASS ] invariant_mechanism_rejection
[ PASS ] invariant_input_validation
----
3/4 tests passed
Start 4: test_soak
4/4 Test #4: test_soak ........................ Passed 0.01 sec
75% tests passed, 1 tests failed out of 4
Label Time Summary:
hidden = 0.02 sec*proc (4 tests)
Total Test time (real) = 0.02 sec
The following tests FAILED:
3 - test_invariants (Failed)
/tmp/tmp.4LkWH2z2RJ/beamfd/src/solver.cpp: In member function 'beamfd::BeamSystem beamfd::BeamSolver::assemble_system() const':
/tmp/tmp.4LkWH2z2RJ/beamfd/src/solver.cpp:83:10: warning: variable 'add_curvature' set but not used [-Wunused-but-set-variable]
83 | auto add_curvature = [&](std::size_t r, std::size_t j, double coef) {
| ^~~~~~~~~~~~~
Errors while running CTest
FAIL: one or more hidden tests failedReproduce this trial: git checkout 2f94510 && PYTHONPATH=src python3 scripts/build_site.py , then open trial/trial_96087325142d4ce8. Re-running the agent live requires EVAL_PLATFORM_ENABLE_OAUTH_SMOKE=1 and is non-deterministic.
Trial trial_96087325142d4ce8 · verifier authoritative; classifier explanatory.