diff --git a/app-proto/src/assembly.rs b/app-proto/src/assembly.rs index 278a8f9..1684716 100644 --- a/app-proto/src/assembly.rs +++ b/app-proto/src/assembly.rs @@ -33,9 +33,8 @@ pub struct 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 + // last time the assembly was realized + column_index: usize } impl Element { @@ -63,7 +62,7 @@ impl Element { representation: create_signal(representation), constraints: create_signal(BTreeSet::default()), serial: serial, - column_index: None + column_index: 0 } } @@ -120,13 +119,6 @@ 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 { @@ -134,16 +126,7 @@ pub struct Assembly { pub elements: Signal>, pub constraints: Signal>, - // 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 - // + // solution variety tangent space pub tangent: Signal, // indexing @@ -155,7 +138,7 @@ impl Assembly { Assembly { elements: create_signal(Slab::new()), constraints: create_signal(Slab::new()), - tangent: create_signal(ConfigSubspace::zero(0)), + tangent: create_signal(ConfigSubspace::zero()), elements_by_id: create_signal(FxHashMap::default()) } } @@ -169,6 +152,13 @@ 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 { @@ -219,7 +209,7 @@ impl Assembly { // index the elements self.elements.update_silent(|elts| { for (index, (_, elt)) in elts.into_iter().enumerate() { - elt.column_index = Some(index); + elt.column_index = index; } }); @@ -231,8 +221,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.unwrap(); - let col = elts[subjects.1].column_index.unwrap(); + let row = elts[subjects.0].column_index; + let col = elts[subjects.1].column_index; gram_to_be.push_sym(row, col, cst.lorentz_prod.get_untracked()); } } @@ -242,7 +232,7 @@ impl Assembly { // Gram matrix let mut guess_to_be = DMatrix::::zeros(5, elts.len()); for (_, elt) in elts { - let index = elt.column_index.unwrap(); + let index = elt.column_index; gram_to_be.push_sym(index, index, 1.0); guess_to_be.set_column(index, &elt.representation.get_clone_untracked()); } @@ -289,7 +279,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.unwrap())) + |rep| rep.set_column(0, &config.column(elt.column_index)) ); } @@ -300,67 +290,26 @@ impl Assembly { // --- 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) { + pub fn deform(&self, element_motions: Vec<(ElementKey, DVectorView)>) { /* 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) { + // 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) { console::log_1(&JsValue::from("The assembly is rigid")); + return; } - // 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 - }); + const ELEMENT_DIM: usize = 5; + let assembly_dim = self.elements.with(|elts| elts.len()); + let mut motion_proj = DMatrix::zeros(ELEMENT_DIM, assembly_dim); // 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; - } + // 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)); } // step each element along the mass shell geodesic that matches its @@ -370,24 +319,13 @@ impl Assembly { // 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) - )) - } - }; + 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)); }); } - // bring the configuration back onto the solution variety. this also - // gets the elements' column indices and the saved tangent space back in - // sync + // bring the configuration back onto the solution variety self.realize(); } } \ No newline at end of file diff --git a/app-proto/src/display.rs b/app-proto/src/display.rs index d8ee23c..673985c 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, ElementMotion}}; +use crate::{AppState, assembly::ElementKey}; fn compile_shader( context: &WebGl2RenderingContext, @@ -341,14 +341,7 @@ pub fn Display() -> View { ]) }; let elt_motion: DVector = time_step * TRANSLATION_SPEED * vel_field * rep; - assembly_for_raf.deform( - vec![ - ElementMotion { - key: sel, - velocity: elt_motion.as_view() - } - ] - ); + assembly_for_raf.deform(vec![(sel, elt_motion.as_view())]); scene_changed.set(true); } } diff --git a/app-proto/src/engine.rs b/app-proto/src/engine.rs index e96ea52..2f06035 100644 --- a/app-proto/src/engine.rs +++ b/app-proto/src/engine.rs @@ -88,17 +88,11 @@ impl PartialMatrix { // --- configuration subspaces --- #[derive(Clone)] -pub struct ConfigSubspace { - assembly_dim: usize, - basis: Vec> -} +pub struct ConfigSubspace(Vec>); impl ConfigSubspace { - pub fn zero(assembly_dim: usize) -> ConfigSubspace { - ConfigSubspace { - assembly_dim: assembly_dim, - basis: Vec::new() - } + pub fn zero() -> ConfigSubspace { + ConfigSubspace(Vec::new()) } // approximate the kernel of a symmetric endomorphism of the configuration @@ -125,31 +119,24 @@ impl ConfigSubspace { format!("Eigenvalues used to find kernel: {}", eig.eigenvalues) )); - ConfigSubspace { - assembly_dim: assembly_dim, - basis: basis.collect() - } + ConfigSubspace(basis.collect()) } pub fn dim(&self) -> usize { - self.basis.len() - } - - pub fn assembly_dim(&self) -> usize { - self.assembly_dim + let ConfigSubspace(basis) = self; + basis.len() } // 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 { - 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() - } + let ConfigSubspace(basis) = self; + basis.into_iter().map( + |b| b.column(column_index).dot(&v) * b + ).sum() } } @@ -338,14 +325,14 @@ pub fn realize_gram( state = better_state; history.backoff_steps.push(backoff_steps); }, - None => return (state.config, ConfigSubspace::zero(assembly_dim), false, history) + None => return (state.config, ConfigSubspace::zero(), false, history) }; } let success = state.loss < tol; let tangent = if success { ConfigSubspace::symmetric_kernel(hess, assembly_dim) } else { - ConfigSubspace::zero(assembly_dim) + ConfigSubspace::zero() }; (state.config, tangent, success, history) }