Close kinematic loops: driven joints move four-bars and sliders (#1141)
Driven joints now solve closed kinematic loops, so a four-bar or slider-crank moves as a whole when one joint is driven.
This commit is contained in:
Mateusz Sych authored and GitHub committed
1 parent
2e6951135a
commit
94faf9ce1b
7 files changed
+1025
-64
No files matched your search
+520
-46
@@ -7,14 +7,17 @@
|
||||
//! `T_B · F_b = T_A · F_a · M(q) · Flip`, so `T_B = T_A · F_a · M(q) · Flip · F_b⁻¹`, where `T`
|
||||
//! are occurrence world transforms, `F` the joint frames in each component's own frame, and `q`
|
||||
//! the joint's values (angles in radians, distances in mm). The solver walks the graph from
|
||||
//! grounded occurrences (and the root, which never moves); joints that close a loop are checked
|
||||
//! and reported when they disagree.
|
||||
//! grounded occurrences (and the root, which never moves). Joints that close a loop (a four-bar
|
||||
//! linkage, a slider-crank) are met by moving the values of the other joints around the loop;
|
||||
//! those that still disagree are reported.
|
||||
|
||||
use std::collections::{BTreeMap, BTreeSet, VecDeque};
|
||||
|
||||
use serde::{Deserialize, Serialize};
|
||||
use solidcraft_geom::Vec3;
|
||||
|
||||
use solidcraft_sketch::linalg;
|
||||
|
||||
use crate::{DocError, Document, IDENTITY, Mat, ModelState, Result, apply_point, apply_vector, mat_inverse, mat_mul};
|
||||
|
||||
/// Joints the solver (and the URDF export) reads at most.
|
||||
@@ -320,6 +323,9 @@ pub struct Solution {
|
||||
/// New world transforms of occurrences the joints place.
|
||||
#[serde(skip)]
|
||||
pub world: BTreeMap<u64, Mat>,
|
||||
/// New values of the joints that moved to close loops (joint id → values).
|
||||
#[serde(skip)]
|
||||
pub values: BTreeMap<u64, Vec<f64>>,
|
||||
/// Joints that disagree with the placement found through others: (joint id, name, why).
|
||||
pub conflicts: Vec<(u64, String, String)>,
|
||||
/// Joints that could not be applied (bad origin…): (joint id, name, why).
|
||||
@@ -330,21 +336,42 @@ pub struct Solution {
|
||||
pub fixed: Vec<u64>,
|
||||
}
|
||||
|
||||
/// Joint values by joint id.
|
||||
type Values = BTreeMap<u64, Vec<f64>>;
|
||||
|
||||
/// How far two placements differ: (translation mm, rotation rad).
|
||||
fn mismatch(a: &Mat, b: &Mat) -> (f64, f64) {
|
||||
let t = Vec3::new(a[3][0] - b[3][0], a[3][1] - b[3][1], a[3][2] - b[3][2]).len();
|
||||
let mut tr = 0.0;
|
||||
for i in 0..3 {
|
||||
for k in 0..3 {
|
||||
tr += a[i][k] * b[i][k];
|
||||
// The turn between them from both its cosine and its sine: the cosine alone loses small
|
||||
// angles to rounding.
|
||||
let c = |i: usize, k: usize| (0..3).map(|n| a[i][n] * b[k][n]).sum::<f64>();
|
||||
let sin2 = Vec3::new(c(1, 2) - c(2, 1), c(2, 0) - c(0, 2), c(0, 1) - c(1, 0)).len();
|
||||
(t, sin2.atan2(c(0, 0) + c(1, 1) + c(2, 2) - 1.0))
|
||||
}
|
||||
|
||||
/// Passes over the motion links: a chain of this many links carries a value through in any order.
|
||||
const MAX_LINK_HOPS: usize = 4;
|
||||
|
||||
/// Motion links in order (a few passes for chains), then the linked values' limits.
|
||||
fn apply_links(asm: &Assembly, vals: &mut Values) {
|
||||
for _ in 0..MAX_LINK_HOPS {
|
||||
for l in &asm.links {
|
||||
let Some(src) = vals.get(&l.a).and_then(|q| q.get(l.ia)).copied() else { continue };
|
||||
if let Some(v) = vals.get_mut(&l.b).and_then(|q| q.get_mut(l.ib)) {
|
||||
*v = l.ratio * src + l.offset;
|
||||
}
|
||||
}
|
||||
}
|
||||
for l in &asm.links {
|
||||
if let (Some(j), Some(q)) = (asm.joints.iter().find(|j| j.id == l.b), vals.get_mut(&l.b)) {
|
||||
clamp(j, q);
|
||||
}
|
||||
}
|
||||
(t, ((tr - 1.0) / 2.0).clamp(-1.0, 1.0).acos())
|
||||
}
|
||||
|
||||
/// Effective values of every joint after motion links and limits.
|
||||
pub fn joint_values(asm: &Assembly) -> BTreeMap<u64, Vec<f64>> {
|
||||
let mut vals: BTreeMap<u64, Vec<f64>> = asm
|
||||
let mut vals: Values = asm
|
||||
.joints
|
||||
.iter()
|
||||
.map(|j| {
|
||||
@@ -353,15 +380,7 @@ pub fn joint_values(asm: &Assembly) -> BTreeMap<u64, Vec<f64>> {
|
||||
(j.id, q)
|
||||
})
|
||||
.collect();
|
||||
// Links in order, a few passes for chains.
|
||||
for _ in 0..4 {
|
||||
for l in &asm.links {
|
||||
let Some(src) = vals.get(&l.a).and_then(|q| q.get(l.ia)).copied() else { continue };
|
||||
if let Some(v) = vals.get_mut(&l.b).and_then(|q| q.get_mut(l.ib)) {
|
||||
*v = l.ratio * src + l.offset;
|
||||
}
|
||||
}
|
||||
}
|
||||
apply_links(asm, &mut vals);
|
||||
for j in &asm.joints {
|
||||
if let Some(q) = vals.get_mut(&j.id) {
|
||||
clamp(j, q);
|
||||
@@ -370,12 +389,450 @@ pub fn joint_values(asm: &Assembly) -> BTreeMap<u64, Vec<f64>> {
|
||||
vals
|
||||
}
|
||||
|
||||
/// The turn and offset from A's joint frame to B's: mated (B's z against A's, as faces meet)
|
||||
/// unless flipped.
|
||||
fn flip_of(j: &Joint) -> Mat {
|
||||
let mut m = rot_z(j.angle);
|
||||
m[3][2] = j.offset;
|
||||
if !j.flip { mat_mul(&m, &[[1.0, 0.0, 0.0, 0.0], [0.0, -1.0, 0.0, 0.0], [0.0, 0.0, -1.0, 0.0], [0.0, 0.0, 0.0, 1.0]]) } else { m }
|
||||
}
|
||||
|
||||
/// B's world transform from A's (or A's from B's, going the other way) for the joint's values.
|
||||
fn through(j: &Joint, (fa, fb): &(Mat, Mat), q: &[f64], from_a: bool, t: &Mat) -> Option<Mat> {
|
||||
let link = mat_mul(&mat_mul(&motion(j.kind, q), &flip_of(j)), &mat_inverse(fb)?);
|
||||
if from_a {
|
||||
Some(mat_mul(&mat_mul(t, fa), &link))
|
||||
} else {
|
||||
// T_A = T_B · link⁻¹ · F_a⁻¹
|
||||
Some(mat_mul(&mat_mul(t, &mat_inverse(&link)?), &mat_inverse(fa)?))
|
||||
}
|
||||
}
|
||||
|
||||
/// The spanning tree the solver places occurrences along, out from the fixed ones.
|
||||
struct Tree<'a> {
|
||||
joints: BTreeMap<u64, &'a Joint>,
|
||||
/// Joint frames in component coordinates (A's, B's).
|
||||
frames: BTreeMap<u64, (Mat, Mat)>,
|
||||
fixed: BTreeMap<u64, Mat>,
|
||||
/// Placed occurrence → (the occurrence it is placed from, the joint, from the joint's A side).
|
||||
parent: BTreeMap<u64, (u64, u64, bool)>,
|
||||
/// Placed occurrences in placement order.
|
||||
order: Vec<u64>,
|
||||
}
|
||||
|
||||
impl Tree<'_> {
|
||||
fn step(&self, joint: u64, from_a: bool, t: &Mat, vals: &Values) -> Option<Mat> {
|
||||
let q = vals.get(&joint).map(Vec::as_slice).unwrap_or(&[]);
|
||||
through(self.joints.get(&joint)?, self.frames.get(&joint)?, q, from_a, t)
|
||||
}
|
||||
|
||||
/// Where every placed occurrence is for these values.
|
||||
fn place(&self, vals: &Values) -> BTreeMap<u64, Mat> {
|
||||
let mut placed = self.fixed.clone();
|
||||
for occ in &self.order {
|
||||
let Some(&(from, joint, from_a)) = self.parent.get(occ) else { continue };
|
||||
if let Some(m) = placed.get(&from).and_then(|t| self.step(joint, from_a, t, vals)) {
|
||||
placed.insert(*occ, m);
|
||||
}
|
||||
}
|
||||
placed
|
||||
}
|
||||
|
||||
/// The fixed occurrence an occurrence is placed from, and the tree joints from it up to that one.
|
||||
fn path(&self, mut occ: u64) -> (u64, Vec<(u64, bool)>) {
|
||||
let mut up = Vec::new();
|
||||
while let Some(&(from, joint, from_a)) = self.parent.get(&occ) {
|
||||
if up.len() > self.order.len() {
|
||||
break;
|
||||
}
|
||||
up.push((joint, from_a));
|
||||
occ = from;
|
||||
}
|
||||
(occ, up)
|
||||
}
|
||||
|
||||
/// One occurrence's world placement for these values.
|
||||
fn world(&self, occ: u64, vals: &Values) -> Option<Mat> {
|
||||
let (root, up) = self.path(occ);
|
||||
let mut t = *self.fixed.get(&root)?;
|
||||
for (joint, from_a) in up.iter().rev() {
|
||||
t = self.step(*joint, *from_a, &t, vals)?;
|
||||
}
|
||||
Some(t)
|
||||
}
|
||||
}
|
||||
|
||||
/// Values the loop solver moves at most in one group of loops, and loop-closing joints it takes
|
||||
/// on at most: larger ones are only checked.
|
||||
const MAX_LOOP_VALUES: usize = 64;
|
||||
const MAX_LOOPS: usize = 100;
|
||||
/// Residual evaluations one solve may spend on all its loops.
|
||||
const MAX_LOOP_WORK: usize = 200_000;
|
||||
|
||||
/// Held joint values by (joint, value index): the value each moves from.
|
||||
pub type Held = BTreeMap<(u64, usize), f64>;
|
||||
|
||||
/// The values that set the `around` joints' values through motion links, followed back through
|
||||
/// every link to them (several links may write one value; `apply_links`' order decides which
|
||||
/// wins on the way) for as many links as `apply_links` passes: all the values reached, and those
|
||||
/// that no link sets (the ones a solve can move).
|
||||
fn link_chains(asm: &Assembly, tree: &Tree, around: &BTreeSet<u64>, linked: &BTreeSet<(u64, usize)>) -> (BTreeSet<(u64, usize)>, Vec<(u64, usize)>) {
|
||||
let (mut feeds, mut starts) = (BTreeSet::new(), Vec::new());
|
||||
let mut frontier: BTreeSet<(u64, usize)> =
|
||||
around.iter().filter_map(|id| tree.joints.get(id)).flat_map(|j| (0..j.kind.dofs().len()).map(|i| (j.id, i))).collect();
|
||||
let mut seen = frontier.clone();
|
||||
for _ in 0..MAX_LINK_HOPS {
|
||||
let next: BTreeSet<(u64, usize)> =
|
||||
asm.links.iter().filter(|l| frontier.contains(&(l.b, l.ib))).map(|l| (l.a, l.ia)).filter(|v| !seen.contains(v)).collect();
|
||||
if next.is_empty() {
|
||||
break;
|
||||
}
|
||||
for &(id, i) in &next {
|
||||
feeds.insert((id, i));
|
||||
if !linked.contains(&(id, i)) && tree.joints.get(&id).is_some_and(|j| i < j.kind.dofs().len()) {
|
||||
starts.push((id, i));
|
||||
}
|
||||
}
|
||||
seen.extend(next.iter().copied());
|
||||
frontier = next;
|
||||
}
|
||||
(feeds, starts)
|
||||
}
|
||||
|
||||
/// The held values among `touches` that move: (joint, index, from, to), `from` within the
|
||||
/// joint's limits. A suppressed joint still sets values through its motion links.
|
||||
fn held_moves(asm: &Assembly, held: &Held, touches: &BTreeSet<(u64, usize)>, vals: &Values) -> Vec<(u64, usize, f64, f64)> {
|
||||
let mut moves = Vec::new();
|
||||
for ((id, i), from) in held.iter().filter(|(v, _)| touches.contains(v)) {
|
||||
let (Some(j), Some(to)) = (asm.joints.iter().find(|j| j.id == *id), vals.get(id).and_then(|q| q.get(*i))) else { continue };
|
||||
let from = match j.limits.get(*i).copied().flatten() {
|
||||
Some((lo, hi)) => from.max(lo.min(hi)).min(hi.max(lo)),
|
||||
None => *from,
|
||||
};
|
||||
if from != *to {
|
||||
moves.push((*id, *i, from, *to));
|
||||
}
|
||||
}
|
||||
moves
|
||||
}
|
||||
|
||||
/// Close the loops, each group of loops that share joints on its own: see `close_group`.
|
||||
/// Returns the moved joints' new stored values (and updates `vals`), and the closing joints of
|
||||
/// the groups whose held values moved but which can't stay closed on the way (nothing in those
|
||||
/// moves; a group that can't close without moving is left to the conflict check).
|
||||
fn close_loops(tree: &Tree, closing: &[u64], asm: &Assembly, held: &Held, vals: &mut Values) -> (Values, BTreeSet<u64>) {
|
||||
let (mut changed, mut stuck) = (BTreeMap::new(), BTreeSet::new());
|
||||
if closing.len() > MAX_LOOPS {
|
||||
return (Values::new(), stuck);
|
||||
}
|
||||
let linked: BTreeSet<(u64, usize)> = asm.links.iter().map(|l| (l.b, l.ib)).collect();
|
||||
// Each loop: its closing joint and the joints on the tree paths between its sides; loops that
|
||||
// share a joint (or a value feeding them through links) are solved together.
|
||||
let mut groups: Vec<(Vec<u64>, BTreeSet<u64>, BTreeSet<(u64, usize)>)> = Vec::new();
|
||||
for id in closing {
|
||||
let Some(j) = tree.joints.get(id) else { continue };
|
||||
// A joint between an occurrence and itself (only a hand-edited file has one) is only checked.
|
||||
if j.a.occurrence == j.b.occurrence {
|
||||
continue;
|
||||
}
|
||||
let side = |occ: u64| -> BTreeSet<u64> { tree.path(occ).1.iter().map(|(k, _)| *k).collect() };
|
||||
let mut around: BTreeSet<u64> = side(j.a.occurrence).symmetric_difference(&side(j.b.occurrence)).copied().collect();
|
||||
around.insert(*id);
|
||||
// The values the loop depends on: its joints' and those its links read. Loops tie together
|
||||
// through a shared value only (a rigid joint on both paths has none; two values of one
|
||||
// joint can feed two independent loops).
|
||||
let mut touches = link_chains(asm, tree, &around, &linked).0;
|
||||
touches.extend(around.iter().filter_map(|id| tree.joints.get(id)).flat_map(|j| (0..j.kind.dofs().len()).map(|i| (j.id, i))));
|
||||
let mut members = vec![*id];
|
||||
groups.retain(|(c, a, t)| {
|
||||
if t.is_disjoint(&touches) {
|
||||
return true;
|
||||
}
|
||||
members.extend(c.iter().copied());
|
||||
around.extend(a.iter().copied());
|
||||
touches.extend(t.iter().copied());
|
||||
false
|
||||
});
|
||||
groups.push((members, around, touches));
|
||||
}
|
||||
let spent = std::cell::Cell::new(0usize);
|
||||
for (group, around, touches) in &groups {
|
||||
let moves = held_moves(asm, held, touches, vals);
|
||||
match close_group(tree, group, around, &moves, asm, held, &linked, &spent, vals) {
|
||||
Some(m) => changed.extend(m),
|
||||
None if !moves.is_empty() => stuck.extend(group.iter().copied()),
|
||||
None => {}
|
||||
}
|
||||
}
|
||||
// The stored values with what each group moved (groups may move different values of one
|
||||
// joint).
|
||||
let mut moved = Values::new();
|
||||
for ((id, i), v) in changed {
|
||||
let Some(j) = tree.joints.get(&id) else { continue };
|
||||
let q = moved.entry(id).or_insert_with(|| {
|
||||
let mut q = j.values.clone();
|
||||
q.resize(j.kind.dofs().len(), 0.0);
|
||||
q
|
||||
});
|
||||
if let Some(slot) = q.get_mut(i) {
|
||||
*slot = v;
|
||||
}
|
||||
}
|
||||
(moved, stuck)
|
||||
}
|
||||
|
||||
/// Close one group of loops: move the values of the joints `around` them (all but the `held`
|
||||
/// values, motion-link targets and rigid joints, within their limits; a linked value moves
|
||||
/// through the start of its link chain) until every loop-closing joint meets. A loop open at the
|
||||
/// start is closed first, from where it is. Then the held values `moves` go from the values they
|
||||
/// are held from to their current ones in small steps, each solved from the last, so the mechanism
|
||||
/// stays on its branch (an open four-bar does not flip crossed). Returns the values that changed
|
||||
/// from what is stored, by (joint, index), and updates `vals` (empty when there is nothing to do
|
||||
/// or the group is too big to solve), or `None` when the loops can't close, or can't stay closed
|
||||
/// on the way.
|
||||
fn close_group(
|
||||
tree: &Tree,
|
||||
closing: &[u64],
|
||||
around: &BTreeSet<u64>,
|
||||
moves: &[(u64, usize, f64, f64)],
|
||||
asm: &Assembly,
|
||||
held: &Held,
|
||||
linked: &BTreeSet<(u64, usize)>,
|
||||
spent: &std::cell::Cell<usize>,
|
||||
vals: &mut Values,
|
||||
) -> Option<BTreeMap<(u64, usize), f64>> {
|
||||
let starts = link_chains(asm, tree, around, linked).1;
|
||||
let mut vars: Vec<(u64, usize)> = Vec::new();
|
||||
for j in around.iter().filter_map(|id| tree.joints.get(id)) {
|
||||
vars.extend((0..j.kind.dofs().len()).map(|i| (j.id, i)).filter(|v| !linked.contains(v) && !held.contains_key(v)));
|
||||
}
|
||||
for v in starts {
|
||||
if !held.contains_key(&v) && !vars.contains(&v) {
|
||||
vars.push(v);
|
||||
}
|
||||
}
|
||||
if vars.len() > MAX_LOOP_VALUES {
|
||||
return Some(BTreeMap::new());
|
||||
}
|
||||
let bounds: Vec<(f64, f64)> = vars
|
||||
.iter()
|
||||
.map(|(id, i)| match tree.joints.get(id).and_then(|j| j.limits.get(*i).copied().flatten()) {
|
||||
Some((lo, hi)) => (lo.min(hi), hi.max(lo)),
|
||||
None => (f64::NEG_INFINITY, f64::INFINITY),
|
||||
})
|
||||
.collect();
|
||||
// Turns weigh like moves through points a loop's size apart.
|
||||
let (mut lo, mut hi) = (Vec3::new(f64::MAX, f64::MAX, f64::MAX), Vec3::new(f64::MIN, f64::MIN, f64::MIN));
|
||||
for id in around {
|
||||
let (Some(j), Some((fa, _))) = (tree.joints.get(id), tree.frames.get(id)) else { continue };
|
||||
let Some(t) = tree.world(j.a.occurrence, vals) else { continue };
|
||||
let p = apply_point(&mat_mul(&t, fa), Vec3::ZERO);
|
||||
lo = Vec3::new(lo.x.min(p.x), lo.y.min(p.y), lo.z.min(p.z));
|
||||
hi = Vec3::new(hi.x.max(p.x), hi.y.max(p.y), hi.z.max(p.z));
|
||||
}
|
||||
let scale = if lo.x <= hi.x { (hi - lo).len().clamp(1.0, 1e6) } else { 1.0 };
|
||||
// Closer than the conflict check needs: 1e-7 mm, and turns of 1e-10 rad across the loop.
|
||||
let tol = (1e-10 * scale).min(1e-7);
|
||||
// The values with the held ones `t` of the way along (and the unknowns at `x`).
|
||||
let place = |work: &mut Values, x: &[f64], t: f64| {
|
||||
for (id, i, from, to) in moves {
|
||||
if let Some(slot) = work.get_mut(id).and_then(|q| q.get_mut(*i)) {
|
||||
*slot = from + (to - from) * t;
|
||||
}
|
||||
}
|
||||
for ((id, i), (v, (lo, hi))) in vars.iter().zip(x.iter().zip(&bounds)) {
|
||||
if let Some(slot) = work.get_mut(id).and_then(|q| q.get_mut(*i)) {
|
||||
*slot = v.max(*lo).min(*hi);
|
||||
}
|
||||
}
|
||||
apply_links(asm, work);
|
||||
};
|
||||
let x0: Vec<f64> = vars.iter().map(|(id, i)| vals.get(id).and_then(|q| q.get(*i)).copied().unwrap_or(0.0)).collect();
|
||||
// Steps that move no value around the loops (links included) more than 5 degrees or a
|
||||
// twentieth of the loops' size.
|
||||
// How many such steps apart two sets of values are, over the loops' values.
|
||||
let stride = |a: &Values, b: &Values| -> f64 {
|
||||
let mut most = 0.0f64;
|
||||
for j in around.iter().filter_map(|id| tree.joints.get(id)) {
|
||||
for i in 0..j.kind.dofs().len() {
|
||||
let (Some(p), Some(q)) = (a.get(&j.id).and_then(|v| v.get(i)), b.get(&j.id).and_then(|v| v.get(i))) else { continue };
|
||||
let step = if j.kind.is_angle(i) { 5f64.to_radians() } else { 0.05 * scale };
|
||||
most = most.max((q - p).abs() / step);
|
||||
}
|
||||
}
|
||||
most
|
||||
};
|
||||
let effective = |x: &[f64], t: f64| {
|
||||
let mut v = vals.clone();
|
||||
place(&mut v, x, t);
|
||||
v
|
||||
};
|
||||
let mut first = 1.0f64;
|
||||
if !moves.is_empty() {
|
||||
first = first.min(1.0 / stride(&effective(&x0, 0.0), &effective(&x0, 1.0)));
|
||||
}
|
||||
let mut work = vals.clone();
|
||||
// How far each closing joint is from meeting with the held values `t` of the way along:
|
||||
// three points of B's joint frame, where the joint wants them against where the tree puts them.
|
||||
let mut gaps = |x: &[f64], t: f64| -> Option<Vec<f64>> {
|
||||
spent.set(spent.get() + 1);
|
||||
if spent.get() > MAX_LOOP_WORK {
|
||||
return None;
|
||||
}
|
||||
place(&mut work, x, t);
|
||||
let mut out = Vec::with_capacity(closing.len() * 9);
|
||||
for id in closing {
|
||||
let (j, (_, fb)) = (tree.joints.get(id)?, tree.frames.get(id)?);
|
||||
let tb = tree.world(j.b.occurrence, &work)?;
|
||||
let want = tree.step(*id, true, &tree.world(j.a.occurrence, &work)?, &work)?;
|
||||
for p in [Vec3::ZERO, Vec3::X * scale, Vec3::Y * scale] {
|
||||
let at = apply_point(fb, p);
|
||||
let d = apply_point(&want, at) - apply_point(&tb, at);
|
||||
out.extend([d.x, d.y, d.z]);
|
||||
}
|
||||
}
|
||||
out.iter().all(|v| v.is_finite()).then_some(out)
|
||||
};
|
||||
// Already closed, with no held value moving. (A moving one is followed all the way: a
|
||||
// rocker driven a full turn ends where it started, but can't get there.)
|
||||
if moves.is_empty() && gaps(&x0, 1.0).is_some_and(|r| largest(&r) <= tol) {
|
||||
return Some(BTreeMap::new());
|
||||
}
|
||||
// A loop open at the start closes first.
|
||||
let mut x = x0;
|
||||
if largest(&gaps(&x, 0.0)?) > tol && !least_squares(&mut |x| gaps(x, 0.0), &mut x, &bounds, tol) {
|
||||
return None;
|
||||
}
|
||||
let (mut t, mut dt) = (if moves.is_empty() { 1.0 } else { 0.0 }, first);
|
||||
let mut at = effective(&x, t);
|
||||
while t < 1.0 {
|
||||
if dt < first * 1e-4 {
|
||||
return None;
|
||||
}
|
||||
let t1 = (t + dt).min(1.0);
|
||||
let mut trial = x.clone();
|
||||
if least_squares(&mut |x| gaps(x, t1), &mut trial, &bounds, tol) {
|
||||
// A step that moved some loop value too far (through a link and a limit, or close to
|
||||
// a dead point) is taken again in halves.
|
||||
let next = effective(&trial, t1);
|
||||
if stride(&at, &next) > 2.0 {
|
||||
dt *= 0.5;
|
||||
continue;
|
||||
}
|
||||
(x, t, dt, at) = (trial, t1, (dt * 2.0).min(first), next);
|
||||
} else if spent.get() > MAX_LOOP_WORK {
|
||||
return None;
|
||||
} else {
|
||||
dt *= 0.5;
|
||||
}
|
||||
}
|
||||
if largest(&gaps(&x, 1.0)?) > tol {
|
||||
return None;
|
||||
}
|
||||
// The values that changed from what is stored.
|
||||
let moved = vars
|
||||
.iter()
|
||||
.zip(&x)
|
||||
.filter(|((id, i), v)| tree.joints.get(id).is_some_and(|j| j.values.get(*i).is_none_or(|s| (s - *v).abs() > 1e-12)))
|
||||
.map(|(var, v)| (*var, *v))
|
||||
.collect();
|
||||
*vals = work;
|
||||
Some(moved)
|
||||
}
|
||||
|
||||
fn largest(r: &[f64]) -> f64 {
|
||||
r.iter().fold(0.0f64, |m, v| m.max(v.abs()))
|
||||
}
|
||||
|
||||
/// Levenberg–Marquardt: move `x` (kept within `bounds`) until every value of `f` is within
|
||||
/// `tol` of zero. False when it can't get there.
|
||||
fn least_squares(f: &mut dyn FnMut(&[f64]) -> Option<Vec<f64>>, x: &mut Vec<f64>, bounds: &[(f64, f64)], tol: f64) -> bool {
|
||||
const H: f64 = 1e-6;
|
||||
let n = x.len();
|
||||
let norm = |r: &[f64]| r.iter().map(|v| v * v).sum::<f64>();
|
||||
let Some(mut r) = f(x) else { return false };
|
||||
if n == 0 {
|
||||
return largest(&r) <= tol;
|
||||
}
|
||||
let mut lambda = 1e-3;
|
||||
for _ in 0..50 {
|
||||
if largest(&r) <= tol {
|
||||
return true;
|
||||
}
|
||||
// The Jacobian by central differences, then the normal equations JᵀJ δ = −Jᵀr.
|
||||
let m = r.len();
|
||||
let mut jac = vec![0.0; m * n];
|
||||
for k in 0..n {
|
||||
let (mut xp, mut xm) = (x.clone(), x.clone());
|
||||
if let (Some(p), Some(q)) = (xp.get_mut(k), xm.get_mut(k)) {
|
||||
*p += H;
|
||||
*q -= H;
|
||||
}
|
||||
let (Some(rp), Some(rm)) = (f(&xp), f(&xm)) else { return false };
|
||||
for (i, (a, b)) in rp.iter().zip(&rm).enumerate() {
|
||||
if let Some(slot) = jac.get_mut(i * n + k) {
|
||||
*slot = (a - b) / (2.0 * H);
|
||||
}
|
||||
}
|
||||
}
|
||||
let mut jtj = linalg::Mat::zeros(n, n);
|
||||
let mut g = vec![0.0; n];
|
||||
for (i, ri) in r.iter().enumerate() {
|
||||
let row = jac.get(i * n..(i + 1) * n).unwrap_or(&[]);
|
||||
for (p, a) in row.iter().enumerate().filter(|(_, a)| **a != 0.0) {
|
||||
if let Some(gp) = g.get_mut(p) {
|
||||
*gp -= a * ri;
|
||||
}
|
||||
for (q, b) in row.iter().enumerate().skip(p) {
|
||||
jtj.add(p, q, a * b);
|
||||
}
|
||||
}
|
||||
}
|
||||
for p in 0..n {
|
||||
for q in 0..p {
|
||||
jtj.set(p, q, jtj.get(q, p));
|
||||
}
|
||||
}
|
||||
let mut stepped = false;
|
||||
while lambda < 1e12 {
|
||||
let mut a = jtj.clone();
|
||||
for p in 0..n {
|
||||
a.add(p, p, lambda * (1.0 + jtj.get(p, p)));
|
||||
}
|
||||
if let Some(d) = linalg::solve(a, g.clone()) {
|
||||
let next: Vec<f64> = x.iter().zip(&d).zip(bounds).map(|((v, dv), (lo, hi))| (v + dv).max(*lo).min(*hi)).collect();
|
||||
if let Some(rn) = f(&next)
|
||||
&& norm(&rn) < norm(&r)
|
||||
{
|
||||
*x = next;
|
||||
r = rn;
|
||||
lambda = (lambda / 3.0).max(1e-9);
|
||||
stepped = true;
|
||||
break;
|
||||
}
|
||||
}
|
||||
lambda *= 4.0;
|
||||
}
|
||||
if !stepped {
|
||||
return false;
|
||||
}
|
||||
}
|
||||
largest(&r) <= tol
|
||||
}
|
||||
|
||||
/// Place occurrences from the joint graph.
|
||||
pub fn solve(doc: &Document, st: &ModelState) -> Solution {
|
||||
solve_holding(doc, st, &BTreeMap::new())
|
||||
}
|
||||
|
||||
/// Place occurrences from the joint graph. Joints that close loops are met by moving the values
|
||||
/// of the other joints around the loop. The `held` values (the ones being driven) keep theirs,
|
||||
/// moving there from the values given (where the loops last closed) so the others follow them.
|
||||
pub fn solve_holding(doc: &Document, st: &ModelState, held: &Held) -> Solution {
|
||||
let asm = &doc.assembly;
|
||||
let mut sol = Solution::default();
|
||||
let joints: Vec<&Joint> = asm.joints.iter().filter(|j| !j.suppressed).take(MAX_JOINTS).collect();
|
||||
let vals = joint_values(asm);
|
||||
let mut vals = joint_values(asm);
|
||||
// Joint frames in component coordinates (errors leave the joint out).
|
||||
let mut frames: BTreeMap<u64, (Mat, Mat)> = BTreeMap::new();
|
||||
for j in &joints {
|
||||
@@ -403,30 +860,15 @@ pub fn solve(doc: &Document, st: &ModelState) -> Solution {
|
||||
dof.insert(j.a.occurrence, 0);
|
||||
}
|
||||
sol.fixed = placed.keys().copied().collect();
|
||||
let flip_m = |j: &Joint| {
|
||||
let mut m = rot_z(j.angle);
|
||||
m[3][2] = j.offset;
|
||||
// Origins mate (B's z against A's, as faces meet) unless flipped.
|
||||
if !j.flip { mat_mul(&m, &[[1.0, 0.0, 0.0, 0.0], [0.0, -1.0, 0.0, 0.0], [0.0, 0.0, -1.0, 0.0], [0.0, 0.0, 0.0, 1.0]]) } else { m }
|
||||
};
|
||||
// B's world transform from A's (or A's from B's, going the other way).
|
||||
let through = |j: &Joint, from_a: bool, t: &Mat| -> Option<Mat> {
|
||||
let (fa, fb) = frames.get(&j.id)?;
|
||||
let q = vals.get(&j.id).cloned().unwrap_or_default();
|
||||
let link = mat_mul(&mat_mul(&motion(j.kind, &q), &flip_m(j)), &mat_inverse(fb)?);
|
||||
if from_a {
|
||||
Some(mat_mul(&mat_mul(t, fa), &link))
|
||||
} else {
|
||||
// T_A = T_B · link⁻¹ · F_a⁻¹
|
||||
Some(mat_mul(&mat_mul(t, &mat_inverse(&link)?), &mat_inverse(fa)?))
|
||||
}
|
||||
};
|
||||
let mut tree =
|
||||
// Every joint's values are known, suppressed ones too (their motion links still set values).
|
||||
Tree { joints: asm.joints.iter().map(|j| (j.id, j)).collect(), frames, fixed: placed.clone(), parent: BTreeMap::new(), order: Vec::new() };
|
||||
let mut used: BTreeSet<u64> = BTreeSet::new();
|
||||
let mut queue: VecDeque<u64> = placed.keys().copied().collect();
|
||||
while let Some(occ) = queue.pop_front() {
|
||||
let Some(t) = placed.get(&occ).copied() else { continue };
|
||||
for j in &joints {
|
||||
if used.contains(&j.id) || !frames.contains_key(&j.id) {
|
||||
if used.contains(&j.id) || !tree.frames.contains_key(&j.id) {
|
||||
continue;
|
||||
}
|
||||
let (other, from_a) = if j.a.occurrence == occ {
|
||||
@@ -439,27 +881,45 @@ pub fn solve(doc: &Document, st: &ModelState) -> Solution {
|
||||
if placed.contains_key(&other) {
|
||||
continue;
|
||||
}
|
||||
let Some(m) = through(j, from_a, &t) else { continue };
|
||||
let Some(m) = tree.step(j.id, from_a, &t, &vals) else { continue };
|
||||
used.insert(j.id);
|
||||
placed.insert(other, m);
|
||||
tree.parent.insert(other, (occ, j.id, from_a));
|
||||
tree.order.push(other);
|
||||
dof.insert(other, dof.get(&occ).copied().unwrap_or(0) + j.kind.dofs().len());
|
||||
queue.push_back(other);
|
||||
}
|
||||
}
|
||||
// Joints that close loops: do they agree?
|
||||
for j in &joints {
|
||||
if used.contains(&j.id) || !frames.contains_key(&j.id) {
|
||||
continue;
|
||||
// Joints that close loops: move the loops' values to meet them, then check they agree.
|
||||
let closing: Vec<u64> = joints
|
||||
.iter()
|
||||
.filter(|j| !used.contains(&j.id) && tree.frames.contains_key(&j.id))
|
||||
.filter(|j| placed.contains_key(&j.a.occurrence) && placed.contains_key(&j.b.occurrence))
|
||||
.map(|j| j.id)
|
||||
.collect();
|
||||
let mut stuck = BTreeSet::new();
|
||||
if !closing.is_empty() {
|
||||
(sol.values, stuck) = close_loops(&tree, &closing, asm, held, &mut vals);
|
||||
if !sol.values.is_empty() {
|
||||
placed = tree.place(&vals);
|
||||
}
|
||||
}
|
||||
for id in &closing {
|
||||
let Some(j) = tree.joints.get(id) else { continue };
|
||||
let (Some(ta), Some(tb)) = (placed.get(&j.a.occurrence), placed.get(&j.b.occurrence)) else { continue };
|
||||
let Some(want) = through(j, true, ta) else { continue };
|
||||
let (dt, da) = mismatch(&want, tb);
|
||||
let (Some(want), Some((_, fb))) = (tree.step(*id, true, ta, &vals), tree.frames.get(id)) else { continue };
|
||||
// Measured at the joint, where the loop solver closes it.
|
||||
let (dt, da) = mismatch(&mat_mul(&want, fb), &mat_mul(tb, fb));
|
||||
if dt > 1e-6 || da > 1e-8 {
|
||||
sol.conflicts.push((
|
||||
j.id,
|
||||
j.name.clone(),
|
||||
format!("off by {dt:.4} mm and {:.3}° from the placement the other joints give", da.to_degrees()),
|
||||
));
|
||||
} else if stuck.contains(id) {
|
||||
// It meets there, but the loop would have to come apart on the way (a rocker
|
||||
// driven a full turn ends where it started).
|
||||
sol.conflicts.push((j.id, j.name.clone(), "the loop can't stay closed on the way to these values".into()));
|
||||
}
|
||||
}
|
||||
for o in &doc.occurrences {
|
||||
@@ -472,10 +932,16 @@ pub fn solve(doc: &Document, st: &ModelState) -> Solution {
|
||||
}
|
||||
|
||||
impl Document {
|
||||
/// Apply a solution: occurrence transforms (in their parents) for the occurrences it places.
|
||||
/// Returns true when something moved.
|
||||
/// Apply a solution: occurrence transforms (in their parents) for the occurrences it places,
|
||||
/// and the values of the joints that moved to close loops. Returns true when something moved.
|
||||
pub fn apply_joint_solution(&mut self, sol: &Solution) -> bool {
|
||||
let mut changed = false;
|
||||
for (id, q) in &sol.values {
|
||||
if let Some(j) = self.assembly.joints.iter_mut().find(|j| j.id == *id) {
|
||||
j.values.clone_from(q);
|
||||
changed = true;
|
||||
}
|
||||
}
|
||||
let parents: BTreeMap<u64, u64> = self.occurrences.iter().map(|o| (o.id, o.parent)).collect();
|
||||
for (occ, world) in &sol.world {
|
||||
let parent = parents.get(occ).copied().unwrap_or(0);
|
||||
@@ -536,4 +1002,12 @@ mod tests {
|
||||
let f = frame(Vec3::new(1.0, 2.0, 3.0), Vec3::Z, None).unwrap();
|
||||
assert!(close(&mat_mul(&f, &mat_inverse(&f).unwrap()), &IDENTITY));
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn mismatch_resolves_tiny_turns() {
|
||||
let (_, a) = mismatch(&rot_z(1e-9), &IDENTITY);
|
||||
assert!((a - 1e-9).abs() < 1e-15, "{a}");
|
||||
let (_, a) = mismatch(&rot_z(3.0), &IDENTITY);
|
||||
assert!((a - 3.0).abs() < 1e-12, "{a}");
|
||||
}
|
||||
}
|
||||
Reference in new issue
Block a user