6df0e855cf
Also, correct the check for whether an element had a column index when we started. The previous revision would've gotten the wrong answer for an element without a column index that appeared more than once in the motion.
393 lines
No EOL
16 KiB
Rust
393 lines
No EOL
16 KiB
Rust
use nalgebra::{DMatrix, DVector, DVectorView, Vector3};
|
|
use rustc_hash::FxHashMap;
|
|
use slab::Slab;
|
|
use std::{collections::BTreeSet, sync::atomic::{AtomicU64, Ordering}};
|
|
use sycamore::prelude::*;
|
|
use web_sys::{console, wasm_bindgen::JsValue}; /* DEBUG */
|
|
|
|
use crate::engine::{realize_gram, ConfigSubspace, PartialMatrix, Q};
|
|
|
|
// the types of the keys we use to access an assembly's elements and constraints
|
|
pub type ElementKey = usize;
|
|
pub type ConstraintKey = usize;
|
|
|
|
pub type ElementColor = [f32; 3];
|
|
|
|
/* KLUDGE */
|
|
// we should reconsider this design when we build a system for switching between
|
|
// assemblies. at that point, we might want to switch to hierarchical keys,
|
|
// where each each element has a key that identifies it within its assembly and
|
|
// each assembly has a key that identifies it within the sesssion
|
|
static NEXT_ELEMENT_SERIAL: AtomicU64 = AtomicU64::new(0);
|
|
|
|
#[derive(Clone, PartialEq)]
|
|
pub struct Element {
|
|
pub id: String,
|
|
pub label: String,
|
|
pub color: ElementColor,
|
|
pub representation: Signal<DVector<f64>>,
|
|
pub constraints: Signal<BTreeSet<ConstraintKey>>,
|
|
|
|
// a serial number, assigned by `Element::new`, that uniquely identifies
|
|
// each element
|
|
pub serial: u64,
|
|
|
|
// the configuration matrix column index that was assigned to this element
|
|
// last time the assembly was realized, or `None` if the element has never
|
|
// been through a realization
|
|
column_index: Option<usize>
|
|
}
|
|
|
|
impl Element {
|
|
pub fn new(
|
|
id: String,
|
|
label: String,
|
|
color: ElementColor,
|
|
representation: DVector<f64>
|
|
) -> Element {
|
|
// take the next serial number, panicking if that was the last number we
|
|
// had left. the technique we use to panic on overflow is taken from
|
|
// _Rust Atomics and Locks_, by Mara Bos
|
|
//
|
|
// https://marabos.nl/atomics/atomics.html#example-handle-overflow
|
|
//
|
|
let serial = NEXT_ELEMENT_SERIAL.fetch_update(
|
|
Ordering::SeqCst, Ordering::SeqCst,
|
|
|serial| serial.checked_add(1)
|
|
).expect("Out of serial numbers for elements");
|
|
|
|
Element {
|
|
id: id,
|
|
label: label,
|
|
color: color,
|
|
representation: create_signal(representation),
|
|
constraints: create_signal(BTreeSet::default()),
|
|
serial: serial,
|
|
column_index: None
|
|
}
|
|
}
|
|
|
|
// the smallest positive depth, represented as a multiple of `dir`, where
|
|
// the line generated by `dir` hits the element (which is assumed to be a
|
|
// sphere). returns `None` if the line misses the sphere. this function
|
|
// should be kept synchronized with `sphere_cast` in `inversive.frag`, which
|
|
// does essentially the same thing on the GPU side
|
|
pub fn cast(&self, dir: Vector3<f64>, assembly_to_world: &DMatrix<f64>) -> Option<f64> {
|
|
// if `a/b` is less than this threshold, we approximate
|
|
// `a*u^2 + b*u + c` by the linear function `b*u + c`
|
|
const DEG_THRESHOLD: f64 = 1e-9;
|
|
|
|
let rep = self.representation.with_untracked(|rep| assembly_to_world * rep);
|
|
let a = -rep[3] * dir.norm_squared();
|
|
let b = rep.rows_range(..3).dot(&dir);
|
|
let c = -rep[4];
|
|
|
|
let adjust = 4.0*a*c/(b*b);
|
|
if adjust < 1.0 {
|
|
// as long as `b` is non-zero, the linear approximation of
|
|
//
|
|
// a*u^2 + b*u + c
|
|
//
|
|
// at `u = 0` will reach zero at a finite depth `u_lin`. the root of
|
|
// the quadratic adjacent to `u_lin` is stored in `lin_root`. if
|
|
// both roots have the same sign, `lin_root` will be the one closer
|
|
// to `u = 0`
|
|
let square_rect_ratio = 1.0 + (1.0 - adjust).sqrt();
|
|
let lin_root = -(2.0*c)/b / square_rect_ratio;
|
|
if a.abs() > DEG_THRESHOLD * b.abs() {
|
|
if lin_root > 0.0 {
|
|
Some(lin_root)
|
|
} else {
|
|
let other_root = -b/(2.*a) * square_rect_ratio;
|
|
(other_root > 0.0).then_some(other_root)
|
|
}
|
|
} else {
|
|
(lin_root > 0.0).then_some(lin_root)
|
|
}
|
|
} else {
|
|
// the line through `dir` misses the sphere completely
|
|
None
|
|
}
|
|
}
|
|
}
|
|
|
|
#[derive(Clone)]
|
|
pub struct Constraint {
|
|
pub subjects: (ElementKey, ElementKey),
|
|
pub lorentz_prod: Signal<f64>,
|
|
pub lorentz_prod_text: Signal<String>,
|
|
pub lorentz_prod_valid: Signal<bool>,
|
|
pub active: Signal<bool>
|
|
}
|
|
|
|
pub struct ElementMotion<'a> {
|
|
pub key: ElementKey,
|
|
pub velocity: DVectorView<'a, f64>
|
|
}
|
|
|
|
type AssemblyMotion<'a> = Vec<ElementMotion<'a>>;
|
|
|
|
// a complete, view-independent description of an assembly
|
|
#[derive(Clone)]
|
|
pub struct Assembly {
|
|
// elements and constraints
|
|
pub elements: Signal<Slab<Element>>,
|
|
pub constraints: Signal<Slab<Constraint>>,
|
|
|
|
// solution variety tangent space. the basis vectors are stored in
|
|
// configuration matrix format, ordered according to the elements' column
|
|
// indices. when you realize the assembly, every element that's present
|
|
// during realization gets a column index and is reflected in the tangent
|
|
// space. since the methods in this module never assign column indices
|
|
// without later realizing the assembly, we get the following invariant:
|
|
//
|
|
// (1) if an element has a column index, its tangent motions can be found
|
|
// in that column of the tangent space basis matrices
|
|
//
|
|
pub tangent: Signal<ConfigSubspace>,
|
|
|
|
// indexing
|
|
pub elements_by_id: Signal<FxHashMap<String, ElementKey>>
|
|
}
|
|
|
|
impl Assembly {
|
|
pub fn new() -> Assembly {
|
|
Assembly {
|
|
elements: create_signal(Slab::new()),
|
|
constraints: create_signal(Slab::new()),
|
|
tangent: create_signal(ConfigSubspace::zero(0)),
|
|
elements_by_id: create_signal(FxHashMap::default())
|
|
}
|
|
}
|
|
|
|
// --- inserting elements and constraints ---
|
|
|
|
// insert an element into the assembly without checking whether we already
|
|
// have an element with the same identifier. any element that does have the
|
|
// same identifier will get kicked out of the `elements_by_id` index
|
|
fn insert_element_unchecked(&self, elt: Element) {
|
|
let id = elt.id.clone();
|
|
let key = self.elements.update(|elts| elts.insert(elt));
|
|
self.elements_by_id.update(|elts_by_id| elts_by_id.insert(id, key));
|
|
}
|
|
|
|
pub fn try_insert_element(&self, elt: Element) -> bool {
|
|
let can_insert = self.elements_by_id.with_untracked(
|
|
|elts_by_id| !elts_by_id.contains_key(&elt.id)
|
|
);
|
|
if can_insert {
|
|
self.insert_element_unchecked(elt);
|
|
}
|
|
can_insert
|
|
}
|
|
|
|
pub fn insert_new_element(&self) {
|
|
// find the next unused identifier in the default sequence
|
|
let mut id_num = 1;
|
|
let mut id = format!("sphere{}", id_num);
|
|
while self.elements_by_id.with_untracked(
|
|
|elts_by_id| elts_by_id.contains_key(&id)
|
|
) {
|
|
id_num += 1;
|
|
id = format!("sphere{}", id_num);
|
|
}
|
|
|
|
// create and insert a new element
|
|
self.insert_element_unchecked(
|
|
Element::new(
|
|
id,
|
|
format!("Sphere {}", id_num),
|
|
[0.75_f32, 0.75_f32, 0.75_f32],
|
|
DVector::<f64>::from_column_slice(&[0.0, 0.0, 0.0, 0.5, -0.5])
|
|
)
|
|
);
|
|
}
|
|
|
|
pub fn insert_constraint(&self, constraint: Constraint) {
|
|
let subjects = constraint.subjects;
|
|
let key = self.constraints.update(|csts| csts.insert(constraint));
|
|
let subject_constraints = self.elements.with(
|
|
|elts| (elts[subjects.0].constraints, elts[subjects.1].constraints)
|
|
);
|
|
subject_constraints.0.update(|csts| csts.insert(key));
|
|
subject_constraints.1.update(|csts| csts.insert(key));
|
|
}
|
|
|
|
// --- realization ---
|
|
|
|
pub fn realize(&self) {
|
|
// index the elements
|
|
self.elements.update_silent(|elts| {
|
|
for (index, (_, elt)) in elts.into_iter().enumerate() {
|
|
elt.column_index = Some(index);
|
|
}
|
|
});
|
|
|
|
// set up the Gram matrix and the initial configuration matrix
|
|
let (gram, guess) = self.elements.with_untracked(|elts| {
|
|
// set up the off-diagonal part of the Gram matrix
|
|
let mut gram_to_be = PartialMatrix::new();
|
|
self.constraints.with_untracked(|csts| {
|
|
for (_, cst) in csts {
|
|
if cst.active.get_untracked() && cst.lorentz_prod_valid.get_untracked() {
|
|
let subjects = cst.subjects;
|
|
let row = elts[subjects.0].column_index.unwrap();
|
|
let col = elts[subjects.1].column_index.unwrap();
|
|
gram_to_be.push_sym(row, col, cst.lorentz_prod.get_untracked());
|
|
}
|
|
}
|
|
});
|
|
|
|
// set up the initial configuration matrix and the diagonal of the
|
|
// Gram matrix
|
|
let mut guess_to_be = DMatrix::<f64>::zeros(5, elts.len());
|
|
for (_, elt) in elts {
|
|
let index = elt.column_index.unwrap();
|
|
gram_to_be.push_sym(index, index, 1.0);
|
|
guess_to_be.set_column(index, &elt.representation.get_clone_untracked());
|
|
}
|
|
|
|
(gram_to_be, guess_to_be)
|
|
});
|
|
|
|
/* DEBUG */
|
|
// log the Gram matrix
|
|
console::log_1(&JsValue::from("Gram matrix:"));
|
|
gram.log_to_console();
|
|
|
|
/* DEBUG */
|
|
// log the initial configuration matrix
|
|
console::log_1(&JsValue::from("Old configuration:"));
|
|
for j in 0..guess.nrows() {
|
|
let mut row_str = String::new();
|
|
for k in 0..guess.ncols() {
|
|
row_str.push_str(format!(" {:>8.3}", guess[(j, k)]).as_str());
|
|
}
|
|
console::log_1(&JsValue::from(row_str));
|
|
}
|
|
|
|
// look for a configuration with the given Gram matrix
|
|
let (config, tangent, success, history) = realize_gram(
|
|
&gram, guess, &[],
|
|
1.0e-12, 0.5, 0.9, 1.1, 200, 110
|
|
);
|
|
|
|
/* DEBUG */
|
|
// report the outcome of the search
|
|
console::log_1(&JsValue::from(
|
|
if success {
|
|
"Target accuracy achieved!"
|
|
} else {
|
|
"Failed to reach target accuracy"
|
|
}
|
|
));
|
|
console::log_2(&JsValue::from("Steps:"), &JsValue::from(history.scaled_loss.len() - 1));
|
|
console::log_2(&JsValue::from("Loss:"), &JsValue::from(*history.scaled_loss.last().unwrap()));
|
|
console::log_2(&JsValue::from("Tangent dimension:"), &JsValue::from(tangent.dim()));
|
|
|
|
if success {
|
|
// read out the solution
|
|
for (_, elt) in self.elements.get_clone_untracked() {
|
|
elt.representation.update(
|
|
|rep| rep.set_column(0, &config.column(elt.column_index.unwrap()))
|
|
);
|
|
}
|
|
|
|
// save the tangent space
|
|
self.tangent.set_silent(tangent);
|
|
}
|
|
}
|
|
|
|
// --- deformation ---
|
|
|
|
// project the given motion to the tangent space of the solution variety and
|
|
// move the assembly along it. the implementation is based on invariant (1)
|
|
// from above and the following additional invariant:
|
|
//
|
|
// (2) if an element is affected by a constraint, it has a column index
|
|
//
|
|
// we have this invariant because the assembly gets realized each time you
|
|
// add a constraint
|
|
pub fn deform(&self, motion: AssemblyMotion) {
|
|
/* KLUDGE */
|
|
// when the tangent space is zero, deformation won't do anything, but
|
|
// the attempt to deform should be registered in the UI. this console
|
|
// message will do for now
|
|
if self.tangent.with(|tan| tan.dim() <= 0 && tan.assembly_dim() > 0) {
|
|
console::log_1(&JsValue::from("The assembly is rigid"));
|
|
}
|
|
|
|
// give a column index to each moving element that doesn't have one yet.
|
|
// this temporarily breaks invariant (1), but the invariant will be
|
|
// restored when we realize the assembly at the end of the deformation.
|
|
// in the process, we find out how many matrix columns we'll need to
|
|
// hold the deformation
|
|
let realized_dim = self.tangent.with(|tan| tan.assembly_dim());
|
|
let motion_dim = self.elements.update_silent(|elts| {
|
|
let mut next_column_index = realized_dim;
|
|
for elt_motion in motion.iter() {
|
|
let moving_elt = &mut elts[elt_motion.key];
|
|
if moving_elt.column_index.is_none() {
|
|
moving_elt.column_index = Some(next_column_index);
|
|
next_column_index += 1;
|
|
}
|
|
}
|
|
next_column_index
|
|
});
|
|
|
|
// project the element motions onto the tangent space of the solution
|
|
// variety and sum them to get a deformation of the whole assembly. the
|
|
// matrix `motion_proj` that holds the deformation has extra columns for
|
|
// any moving elements that aren't reflected in the saved tangent space
|
|
const ELEMENT_DIM: usize = 5;
|
|
let mut motion_proj = DMatrix::zeros(ELEMENT_DIM, motion_dim);
|
|
for elt_motion in motion {
|
|
// we can unwrap the column index because we know that every moving
|
|
// element has one at this point
|
|
let column_index = self.elements.with_untracked(
|
|
|elts| elts[elt_motion.key].column_index.unwrap()
|
|
);
|
|
|
|
if column_index < realized_dim {
|
|
// this element had a column index when we started, so by
|
|
// invariant (1), it's reflected in the tangent space
|
|
let mut target_columns = motion_proj.columns_mut(0, realized_dim);
|
|
target_columns += self.tangent.with(
|
|
|tan| tan.proj(&elt_motion.velocity, column_index)
|
|
);
|
|
} else {
|
|
// this element didn't have a column index when we started, so
|
|
// by invariant (2), it's unconstrained
|
|
let mut target_column = motion_proj.column_mut(column_index);
|
|
target_column += elt_motion.velocity;
|
|
}
|
|
}
|
|
|
|
// step each element along the mass shell geodesic that matches its
|
|
// velocity in the deformation found above
|
|
/* KLUDGE */
|
|
// since our test assemblies only include spheres, we assume that every
|
|
// element is on the 1 mass shell
|
|
for (_, elt) in self.elements.get_clone_untracked() {
|
|
elt.representation.update_silent(|rep| {
|
|
match elt.column_index {
|
|
Some(column_index) => {
|
|
let rep_next = &*rep + motion_proj.column(column_index);
|
|
let normalizer = rep_next.dot(&(&*Q * &rep_next));
|
|
rep.set_column(0, &(rep_next / normalizer));
|
|
},
|
|
None => {
|
|
console::log_1(&JsValue::from(
|
|
format!("No velocity to unpack for fresh element \"{}\"", elt.id)
|
|
))
|
|
}
|
|
};
|
|
});
|
|
}
|
|
|
|
// bring the configuration back onto the solution variety. this also
|
|
// gets the elements' column indices and the saved tangent space back in
|
|
// sync
|
|
self.realize();
|
|
}
|
|
} |