diff --git a/app-proto/src/assembly.rs b/app-proto/src/assembly.rs index 1684716..278a8f9 100644 --- a/app-proto/src/assembly.rs +++ b/app-proto/src/assembly.rs @@ -33,8 +33,9 @@ pub struct Element { pub serial: u64, // the configuration matrix column index that was assigned to this element - // last time the assembly was realized - column_index: usize + // last time the assembly was realized, or `None` if the element has never + // been through a realization + column_index: Option } impl Element { @@ -62,7 +63,7 @@ impl Element { representation: create_signal(representation), constraints: create_signal(BTreeSet::default()), serial: serial, - column_index: 0 + column_index: None } } @@ -119,6 +120,13 @@ pub struct Constraint { pub active: Signal } +pub struct ElementMotion<'a> { + pub key: ElementKey, + pub velocity: DVectorView<'a, f64> +} + +type AssemblyMotion<'a> = Vec>; + // a complete, view-independent description of an assembly #[derive(Clone)] pub struct Assembly { @@ -126,7 +134,16 @@ pub struct Assembly { pub elements: Signal>, pub constraints: Signal>, - // solution variety tangent space + // 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, // indexing @@ -138,7 +155,7 @@ impl Assembly { Assembly { elements: create_signal(Slab::new()), constraints: create_signal(Slab::new()), - tangent: create_signal(ConfigSubspace::zero()), + tangent: create_signal(ConfigSubspace::zero(0)), elements_by_id: create_signal(FxHashMap::default()) } } @@ -152,13 +169,6 @@ impl Assembly { 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)); - - // realize to update the tangent space - /* KLUDGE */ - // since the newly inserted element is unconstrained, we should be able - // to update the tangent space without recomputing the Hessian and its - // eigendecomposition - self.realize(); } pub fn try_insert_element(&self, elt: Element) -> bool { @@ -209,7 +219,7 @@ impl Assembly { // index the elements self.elements.update_silent(|elts| { for (index, (_, elt)) in elts.into_iter().enumerate() { - elt.column_index = index; + elt.column_index = Some(index); } }); @@ -221,8 +231,8 @@ impl Assembly { 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; - let col = elts[subjects.1].column_index; + 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()); } } @@ -232,7 +242,7 @@ impl Assembly { // Gram matrix let mut guess_to_be = DMatrix::::zeros(5, elts.len()); for (_, elt) in elts { - let index = elt.column_index; + 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()); } @@ -279,7 +289,7 @@ impl Assembly { // 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)) + |rep| rep.set_column(0, &config.column(elt.column_index.unwrap())) ); } @@ -290,26 +300,67 @@ impl Assembly { // --- deformation --- - pub fn deform(&self, element_motions: Vec<(ElementKey, DVectorView)>) { + // 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, we currently need to avoid calling - // its `proj` method, because it will panic rather than returning zero. - // in the future, we'll want a more intentionally designed system for - // handling this case - if self.tangent.with(|tan| tan.dim() <= 0) { + // 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")); - return; } - const ELEMENT_DIM: usize = 5; - let assembly_dim = self.elements.with(|elts| elts.len()); - let mut motion_proj = DMatrix::zeros(ELEMENT_DIM, assembly_dim); + // 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 - for (elt_key, elt_motion) in element_motions { - let column_index = self.elements.with(|elts| elts[elt_key].column_index); - motion_proj += self.tangent.with(|tan| tan.proj(&elt_motion, column_index)); + // 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 @@ -319,13 +370,24 @@ impl Assembly { // element is on the 1 mass shell for (_, elt) in self.elements.get_clone_untracked() { elt.representation.update_silent(|rep| { - let rep_next = &*rep + motion_proj.column(elt.column_index); - let normalizer = rep_next.dot(&(&*Q * &rep_next)); - rep.set_column(0, &(rep_next / normalizer)); + 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 + // 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(); } } \ No newline at end of file diff --git a/app-proto/src/display.rs b/app-proto/src/display.rs index 673985c..d8ee23c 100644 --- a/app-proto/src/display.rs +++ b/app-proto/src/display.rs @@ -14,7 +14,7 @@ use web_sys::{ wasm_bindgen::{JsCast, JsValue} }; -use crate::{AppState, assembly::ElementKey}; +use crate::{AppState, assembly::{ElementKey, ElementMotion}}; fn compile_shader( context: &WebGl2RenderingContext, @@ -341,7 +341,14 @@ pub fn Display() -> View { ]) }; let elt_motion: DVector = time_step * TRANSLATION_SPEED * vel_field * rep; - assembly_for_raf.deform(vec![(sel, elt_motion.as_view())]); + assembly_for_raf.deform( + vec![ + ElementMotion { + key: sel, + velocity: elt_motion.as_view() + } + ] + ); scene_changed.set(true); } } diff --git a/app-proto/src/engine.rs b/app-proto/src/engine.rs index 2f06035..e96ea52 100644 --- a/app-proto/src/engine.rs +++ b/app-proto/src/engine.rs @@ -88,11 +88,17 @@ impl PartialMatrix { // --- configuration subspaces --- #[derive(Clone)] -pub struct ConfigSubspace(Vec>); +pub struct ConfigSubspace { + assembly_dim: usize, + basis: Vec> +} impl ConfigSubspace { - pub fn zero() -> ConfigSubspace { - ConfigSubspace(Vec::new()) + pub fn zero(assembly_dim: usize) -> ConfigSubspace { + ConfigSubspace { + assembly_dim: assembly_dim, + basis: Vec::new() + } } // approximate the kernel of a symmetric endomorphism of the configuration @@ -119,24 +125,31 @@ impl ConfigSubspace { format!("Eigenvalues used to find kernel: {}", eig.eigenvalues) )); - ConfigSubspace(basis.collect()) + ConfigSubspace { + assembly_dim: assembly_dim, + basis: basis.collect() + } } pub fn dim(&self) -> usize { - let ConfigSubspace(basis) = self; - basis.len() + self.basis.len() + } + + pub fn assembly_dim(&self) -> usize { + self.assembly_dim } // find the projection onto this subspace of the motion where the element // with the given column index has velocity `v` - /* TO DO */ - // for the zero subspace, this method's behavior doesn't match its name: it - // panics rather than returning zero pub fn proj(&self, v: &DVectorView, column_index: usize) -> DMatrix { - let ConfigSubspace(basis) = self; - basis.into_iter().map( - |b| b.column(column_index).dot(&v) * b - ).sum() + if self.dim() == 0 { + const ELEMENT_DIM: usize = 5; + DMatrix::zeros(ELEMENT_DIM, self.assembly_dim) + } else { + self.basis.iter().map( + |b| b.column(column_index).dot(&v) * b + ).sum() + } } } @@ -325,14 +338,14 @@ pub fn realize_gram( state = better_state; history.backoff_steps.push(backoff_steps); }, - None => return (state.config, ConfigSubspace::zero(), false, history) + None => return (state.config, ConfigSubspace::zero(assembly_dim), false, history) }; } let success = state.loss < tol; let tangent = if success { ConfigSubspace::symmetric_kernel(hess, assembly_dim) } else { - ConfigSubspace::zero() + ConfigSubspace::zero(assembly_dim) }; (state.config, tangent, success, history) }