Oregami
Repositories/oxedyne/fe2o3

oxedyne/fe2o3/fe2o3_geom/src/rigid.rs

11.7 KiB, 1 run

created by r1870400018:35651, which is this file's identity for as long as the history lasts, whatever it is later renamed to

download · who wrote it · its history

1//! Generic (combinatorial) rigidity of a 2-D bar-joint framework.
2//!
3//! A bar-joint framework is a set of joints in the plane connected by rigid bars. Its
4//! *generic* rigidity -- whether it can flex when the joints sit in general position --
5//! depends only on the underlying graph, not on the exact coordinates. Laman's theorem
6//! makes this precise: a graph on `n` joints is minimally rigid in the plane exactly when
7//! it has `2n - 3` bars and every `k`-joint subset spans at most `2k - 3` bars. The bars
8//! that are independent in this sense form the *2-D rigidity matroid*, and its rank is what
9//! decides rigidity.
10//!
11//! This module computes that rank with the `(2, 3)`-pebble game (Lee and Streinu, "Pebble
12//! game algorithms and sparse graphs"), which tests each bar for independence in turn and
13//! runs in low polynomial time without ever forming the rigidity matrix. From the rank it
14//! reports the internal degrees of freedom and the count of redundant bars.
15//!
16//! The joints are taken as [`Pt`] so a caller passes its geometry directly, but the generic
17//! result reads only the joint *count* and the bar list: two frameworks with the same graph
18//! share a verdict whatever their coordinates. Special positions (three collinear joints,
19//! say) can lose rank in reality, and this generic test deliberately does not model that --
20//! it answers the question a truss puzzle asks, where any two valid layouts are equivalent.
21
22use crate::planar::Pt;
23
24use oxedyne_fe2o3_core::prelude::*;
25
26// The plane gives each joint two pebbles (its two translational freedoms), and a bar is
27// admitted once its endpoints between them hold l + 1 = 4 pebbles, where l = 3 is the count
28// of rigid-body motions of the plane (two translations and one rotation).
29const PEBBLES_PER_JOINT: usize = 2; // k, the spatial dimension
30const ADMIT_THRESHOLD: usize = 4; // l + 1
31
32/// The rigidity verdict for a 2-D bar-joint framework.
33///
34/// `dof` counts the ways the framework can flex beyond the three rigid-body motions of the
35/// plane; it is zero exactly when the framework is rigid. `excess` counts bars that add no
36/// stiffness -- remove any one and the rank is unchanged; it is zero exactly when the
37/// framework carries no redundancy. A framework is minimally rigid when both are zero.
38pub struct Rigidity {
39 pub rank: usize, // independent bars: the rank of the 2-D rigidity matroid
40 pub dof: usize, // internal degrees of freedom, 2n - 3 - rank
41 pub excess: usize, // redundant bars, present bars - rank
42}
43
44impl Rigidity {
45 /// Is the framework rigid (no internal freedom)?
46 pub fn is_rigid(&self) -> bool {
47 self.dof == 0
48 }
49
50 /// Is the framework minimally rigid (rigid and free of redundant bars)?
51 pub fn is_minimally_rigid(&self) -> bool {
52 self.dof == 0 && self.excess == 0
53 }
54
55 /// Returns the combined shortfall, `dof + excess`, which is zero exactly for a minimally
56 /// rigid framework. A single scalar lets a caller gate on one comparison.
57 pub fn flaw(&self) -> usize {
58 self.dof + self.excess
59 }
60}
61
62/// Computes the generic 2-D rigidity of the framework whose joints are `nodes` and whose
63/// bars are the undirected `edges`, each an index pair into `nodes`.
64///
65/// The coordinates carried by `nodes` do not affect the result; only the joint count and the
66/// bar list do (see the module header). Bars may repeat -- a second bar between the same two
67/// joints is always redundant and shows up in `excess`.
68///
69/// Fails when a bar names a joint index outside `nodes`, or joins a joint to itself.
70pub fn analyse(nodes: &[Pt], edges: &[(usize, usize)]) -> Outcome<Rigidity> {
71 let n = nodes.len();
72 for &(a, b) in edges {
73 if a >= n || b >= n {
74 return Err(err!(
75 "Bar ({}, {}) references a joint index outside the {} joints supplied.",
76 a, b, n;
77 Invalid, Input, Range));
78 }
79 if a == b {
80 return Err(err!(
81 "Bar joins joint {} to itself, which is not a bar.", a;
82 Invalid, Input));
83 }
84 }
85
86 let mut game = PebbleGame::new(n);
87 let mut rank = 0;
88 for &(a, b) in edges {
89 if game.admit(a, b) {
90 rank += 1;
91 }
92 }
93
94 let present = edges.len();
95 let budget = (2 * n).saturating_sub(3); // 2n - 3, the minimally rigid bar count
96 let dof = budget.saturating_sub(rank); // rank <= budget for n >= 2, exact there
97 let excess = present.saturating_sub(rank); // rank <= present always, exact
98
99 Ok(Rigidity { rank, dof, excess })
100}
101
102/// The `(2, 3)`-pebble game over a directed orientation of the accepted bars.
103///
104/// The invariant is `pebbles[v] = PEBBLES_PER_JOINT - outdegree(v)`: each joint starts with
105/// its full complement and spends one pebble for each bar oriented out of it. A bar is
106/// independent exactly when its endpoints can gather [`ADMIT_THRESHOLD`] pebbles between
107/// them, pulling free pebbles inwards by reversing directed paths.
108struct PebbleGame {
109 pebbles: Vec<usize>, // free pebbles at each joint
110 out: Vec<Vec<usize>>, // out-neighbours in the current orientation
111}
112
113impl PebbleGame {
114 fn new(n: usize) -> Self {
115 Self {
116 pebbles: vec![PEBBLES_PER_JOINT; n],
117 out: vec![Vec::new(); n],
118 }
119 }
120
121 /// Tries to admit the bar `(u, v)`, returning whether it was independent.
122 ///
123 /// On acceptance the bar is oriented out of an endpoint that holds a pebble, spending it.
124 fn admit(&mut self, u: usize, v: usize) -> bool {
125 while self.pebbles[u] + self.pebbles[v] < ADMIT_THRESHOLD {
126 // Pull one external pebble to either endpoint; stop when neither can find one.
127 if !self.collect(u, u, v) && !self.collect(v, u, v) {
128 break;
129 }
130 }
131 if self.pebbles[u] + self.pebbles[v] >= ADMIT_THRESHOLD {
132 if self.pebbles[u] > 0 {
133 self.pebbles[u] -= 1;
134 self.out[u].push(v);
135 } else {
136 self.pebbles[v] -= 1;
137 self.out[v].push(u);
138 }
139 true
140 } else {
141 false
142 }
143 }
144
145 /// Reverses a directed path from `from` to some free pebble, bringing one pebble to
146 /// `from`. The two endpoints of the bar under test, `ban_a` and `ban_b`, are excluded as
147 /// sources so that each success strictly raises the pebble total on those endpoints.
148 /// Returns whether a pebble was brought.
149 fn collect(&mut self, from: usize, ban_a: usize, ban_b: usize) -> bool {
150 let n = self.pebbles.len();
151 let mut parent = vec![usize::MAX; n];
152 let mut seen = vec![false; n];
153 let mut stack = vec![from];
154 seen[from] = true;
155 let mut found = None;
156 while let Some(v) = stack.pop() {
157 if self.pebbles[v] > 0 && v != ban_a && v != ban_b {
158 found = Some(v);
159 break;
160 }
161 for &w in &self.out[v] {
162 if !seen[w] {
163 seen[w] = true;
164 parent[w] = v;
165 stack.push(w);
166 }
167 }
168 }
169 match found {
170 None => false,
171 Some(target) => {
172 // The pebble moves from `target` to `from`; the path between them flips, so
173 // every intermediate joint keeps its outdegree and only the ends change.
174 self.pebbles[target] -= 1;
175 self.pebbles[from] += 1;
176 let mut w = target;
177 while w != from {
178 let p = parent[w];
179 if let Some(pos) = self.out[p].iter().position(|&z| z == w) {
180 self.out[p].swap_remove(pos);
181 }
182 self.out[w].push(p);
183 w = p;
184 }
185 true
186 },
187 }
188 }
189}
190
191#[cfg(test)]
192mod test {
193 use super::*;
194
195 // A grid of joint positions; the coordinates are immaterial to the generic result, so
196 // any general-position layout serves.
197 fn joints(n: usize) -> Vec<Pt> {
198 (0..n).map(|i| Pt::new(i as f64, (i * i) as f64)).collect()
199 }
200
201 #[test]
202 fn triangle_is_minimally_rigid() -> Outcome<()> {
203 // Three joints, three bars: the smallest rigid framework.
204 let r = res!(analyse(&joints(3), &[(0, 1), (1, 2), (0, 2)]));
205 assert_eq!(r.rank, 3);
206 assert_eq!(r.dof, 0);
207 assert_eq!(r.excess, 0);
208 assert!(r.is_minimally_rigid());
209 Ok(())
210 }
211
212 #[test]
213 fn quadrilateral_is_a_mechanism() -> Outcome<()> {
214 // A four-bar loop has 2n - 4 bars and flexes with one degree of freedom.
215 let r = res!(analyse(&joints(4), &[(0, 1), (1, 2), (2, 3), (3, 0)]));
216 assert_eq!(r.rank, 4);
217 assert_eq!(r.dof, 1);
218 assert_eq!(r.excess, 0);
219 assert!(!r.is_rigid());
220 Ok(())
221 }
222
223 #[test]
224 fn quadrilateral_with_one_diagonal_is_minimally_rigid() -> Outcome<()> {
225 // The single diagonal triangulates the loop: 2n - 3 = 5 bars, no freedom, no waste.
226 let r = res!(analyse(&joints(4), &[(0, 1), (1, 2), (2, 3), (3, 0), (0, 2)]));
227 assert_eq!(r.rank, 5);
228 assert_eq!(r.dof, 0);
229 assert_eq!(r.excess, 0);
230 assert!(r.is_minimally_rigid());
231 Ok(())
232 }
233
234 #[test]
235 fn quadrilateral_with_both_diagonals_is_over_braced() -> Outcome<()> {
236 // The second diagonal adds no stiffness: still rigid, but one redundant bar.
237 let r = res!(analyse(
238 &joints(4),
239 &[(0, 1), (1, 2), (2, 3), (3, 0), (0, 2), (1, 3)],
240 ));
241 assert_eq!(r.rank, 5);
242 assert_eq!(r.dof, 0);
243 assert_eq!(r.excess, 1);
244 assert!(r.is_rigid());
245 assert!(!r.is_minimally_rigid());
246 Ok(())
247 }
248
249 #[test]
250 fn repeated_bar_is_redundant() -> Outcome<()> {
251 // A doubled bar between the same joints can never add stiffness.
252 let r = res!(analyse(&joints(3), &[(0, 1), (1, 2), (0, 2), (0, 1)]));
253 assert_eq!(r.rank, 3);
254 assert_eq!(r.dof, 0);
255 assert_eq!(r.excess, 1);
256 Ok(())
257 }
258
259 #[test]
260 fn right_bar_count_wrong_distribution_is_not_rigid() -> Outcome<()> {
261 // The adversarial case that defeats a bare edge count. Six joints and 2n - 3 = 9
262 // bars, so counting alone would declare it minimally rigid. But the bars are
263 // misplaced: joints 0..3 form a K4 (six bars where five suffice, one redundant),
264 // while the triangle on joints 3, 4, 5 hangs off joint 3 by a single shared joint --
265 // a hinge that lets the two rigid pieces rotate about one another. The framework is a
266 // one-freedom mechanism carrying a redundant bar, not a rigid structure.
267 let edges = [
268 // K4 on {0, 1, 2, 3}: over-braced.
269 (0, 1), (0, 2), (0, 3), (1, 2), (1, 3), (2, 3),
270 // Triangle on {3, 4, 5}: rigid, but joined only at joint 3.
271 (3, 4), (4, 5), (3, 5),
272 ];
273 assert_eq!(edges.len(), 2 * 6 - 3); // the count-only trap: exactly 2n - 3
274 let r = res!(analyse(&joints(6), &edges));
275 assert_eq!(r.rank, 8);
276 assert_eq!(r.dof, 1); // Laman catches the hidden mechanism
277 assert_eq!(r.excess, 1); // and the hidden redundancy
278 assert!(!r.is_rigid());
279 assert!(!r.is_minimally_rigid());
280 Ok(())
281 }
282
283 #[test]
284 fn out_of_range_bar_is_refused() {
285 // Joint index 5 does not exist among three joints.
286 let res = analyse(&joints(3), &[(0, 1), (1, 5)]);
287 assert!(res.is_err());
288 }
289
290 #[test]
291 fn self_bar_is_refused() {
292 // A bar from a joint to itself is not a bar.
293 let res = analyse(&joints(3), &[(0, 0)]);
294 assert!(res.is_err());
295 }
296}