From 2c55a63a6ff959bb389af93e5677815cf4b87117 Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Fri, 6 Dec 2024 14:35:30 -0800 Subject: [PATCH 1/8] Engine: Find the tangent space of the solution variety At the end of the realization routine, use the computed Hessian to find the tangent space of the solution variety, and return it alongside the realization. Since altering the constraints can change the tangent space without changing the solution, we compute the tangent space even when the guess passed to the realization routine is already a solution. --- app-proto/examples/irisawa-hexlet.rs | 2 +- app-proto/examples/point-on-sphere.rs | 2 +- app-proto/examples/three-spheres.rs | 2 +- app-proto/src/assembly.rs | 12 ++- app-proto/src/engine.rs | 131 +++++++++++++++++++++++--- 5 files changed, 129 insertions(+), 20 deletions(-) diff --git a/app-proto/examples/irisawa-hexlet.rs b/app-proto/examples/irisawa-hexlet.rs index fc14f91..2bc94c0 100644 --- a/app-proto/examples/irisawa-hexlet.rs +++ b/app-proto/examples/irisawa-hexlet.rs @@ -2,7 +2,7 @@ use dyna3::engine::{Q, irisawa::realize_irisawa_hexlet}; fn main() { const SCALED_TOL: f64 = 1.0e-12; - let (config, success, history) = realize_irisawa_hexlet(SCALED_TOL); + let (config, _, success, history) = realize_irisawa_hexlet(SCALED_TOL); print!("\nCompleted Gram matrix:{}", config.tr_mul(&*Q) * &config); if success { println!("Target accuracy achieved!"); diff --git a/app-proto/examples/point-on-sphere.rs b/app-proto/examples/point-on-sphere.rs index 0e4765a..13040e5 100644 --- a/app-proto/examples/point-on-sphere.rs +++ b/app-proto/examples/point-on-sphere.rs @@ -18,7 +18,7 @@ fn main() { ]); let frozen = [(3, 0)]; println!(); - let (config, success, history) = realize_gram( + let (config, _, success, history) = realize_gram( &gram, guess, &frozen, 1.0e-12, 0.5, 0.9, 1.1, 200, 110 ); diff --git a/app-proto/examples/three-spheres.rs b/app-proto/examples/three-spheres.rs index d348b18..19acfd1 100644 --- a/app-proto/examples/three-spheres.rs +++ b/app-proto/examples/three-spheres.rs @@ -21,7 +21,7 @@ fn main() { ]) }; println!(); - let (config, success, history) = realize_gram( + let (config, _, success, history) = realize_gram( &gram, guess, &[], 1.0e-12, 0.5, 0.9, 1.1, 200, 110 ); diff --git a/app-proto/src/assembly.rs b/app-proto/src/assembly.rs index fb5bbf7..0d63e45 100644 --- a/app-proto/src/assembly.rs +++ b/app-proto/src/assembly.rs @@ -5,7 +5,7 @@ use std::{collections::BTreeSet, sync::atomic::{AtomicU64, Ordering}}; use sycamore::prelude::*; use web_sys::{console, wasm_bindgen::JsValue}; /* DEBUG */ -use crate::engine::{realize_gram, PartialMatrix}; +use crate::engine::{realize_gram, ConfigSubspace, PartialMatrix}; // the types of the keys we use to access an assembly's elements and constraints pub type ElementKey = usize; @@ -109,7 +109,6 @@ impl Element { } } } - #[derive(Clone)] pub struct Constraint { @@ -127,6 +126,9 @@ pub struct Assembly { pub elements: Signal>, pub constraints: Signal>, + // solution variety tangent space + pub tangent: Signal, + // indexing pub elements_by_id: Signal> } @@ -136,6 +138,7 @@ impl Assembly { Assembly { elements: create_signal(Slab::new()), constraints: create_signal(Slab::new()), + tangent: create_signal(ConfigSubspace::zero()), elements_by_id: create_signal(FxHashMap::default()) } } @@ -247,7 +250,7 @@ impl Assembly { } // look for a configuration with the given Gram matrix - let (config, success, history) = realize_gram( + let (config, tangent, success, history) = realize_gram( &gram, guess, &[], 1.0e-12, 0.5, 0.9, 1.1, 200, 110 ); @@ -271,6 +274,9 @@ impl Assembly { |rep| rep.set_column(0, &config.column(elt.column_index)) ); } + + // save the tangent space + self.tangent.set_silent(tangent); } } } \ No newline at end of file diff --git a/app-proto/src/engine.rs b/app-proto/src/engine.rs index 285a9c6..6a62579 100644 --- a/app-proto/src/engine.rs +++ b/app-proto/src/engine.rs @@ -1,5 +1,5 @@ use lazy_static::lazy_static; -use nalgebra::{Const, DMatrix, DVector, Dyn}; +use nalgebra::{Const, DMatrix, DVector, DVectorView, Dyn, SymmetricEigen}; use web_sys::{console, wasm_bindgen::JsValue}; /* DEBUG */ // --- elements --- @@ -85,6 +85,47 @@ impl PartialMatrix { } } +// --- configuration subspaces --- + +pub struct ConfigSubspace(Vec>); + +impl ConfigSubspace { + pub fn zero() -> ConfigSubspace { + ConfigSubspace(Vec::new()) + } + + // approximate the kernel of a symmetric endomorphism of the configuration + // space for `assembly_dim` elements. we consider an eigenvector to be part + // of the kernel if its eigenvalue is smaller than the constant `THRESHOLD` + fn symmetric_kernel(a: DMatrix, assembly_dim: usize) -> ConfigSubspace { + const ELEMENT_DIM: usize = 5; + const THRESHOLD: f64 = 1.0e-9; + let eig = SymmetricEigen::new(a); + let eig_vecs = eig.eigenvectors.column_iter(); + let eig_pairs = eig.eigenvalues.iter().zip(eig_vecs); + let basis = eig_pairs.filter_map( + |(λ, v)| (λ.abs() < THRESHOLD).then_some( + Into::>::into( + v.reshape_generic(Dyn(ELEMENT_DIM), Dyn(assembly_dim)) + ) + ) + ); + ConfigSubspace(basis.collect()) + } + + // 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 + 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() + } +} + // --- descent history --- pub struct DescentHistory { @@ -181,7 +222,7 @@ pub fn realize_gram( reg_scale: f64, max_descent_steps: i32, max_backoff_steps: i32 -) -> (DMatrix, bool, DescentHistory) { +) -> (DMatrix, ConfigSubspace, bool, DescentHistory) { // start the descent history let mut history = DescentHistory::new(); @@ -201,12 +242,8 @@ pub fn realize_gram( // use Newton's method with backtracking and gradient descent backup let mut state = SearchState::from_config(gram, guess); + let mut hess = DMatrix::zeros(element_dim, assembly_dim); for _ in 0..max_descent_steps { - // stop if the loss is tolerably low - history.config.push(state.config.clone()); - history.scaled_loss.push(state.loss / scale_adjustment); - if state.loss < tol { break; } - // find the negative gradient of the loss function let neg_grad = 4.0 * &*Q * &state.config * &state.err_proj; let mut neg_grad_stacked = neg_grad.clone().reshape_generic(Dyn(total_dim), Const::<1>); @@ -229,7 +266,7 @@ pub fn realize_gram( hess_cols.push(deriv_grad.reshape_generic(Dyn(total_dim), Const::<1>)); } } - let mut hess = DMatrix::from_columns(hess_cols.as_slice()); + hess = DMatrix::from_columns(hess_cols.as_slice()); // regularize the Hessian let min_eigval = hess.symmetric_eigenvalues().min(); @@ -249,6 +286,11 @@ pub fn realize_gram( hess[(k, k)] = 1.0; } + // stop if the loss is tolerably low + history.config.push(state.config.clone()); + history.scaled_loss.push(state.loss / scale_adjustment); + if state.loss < tol { break; } + // compute the Newton step /* we need to either handle or eliminate the case where the minimum @@ -256,7 +298,7 @@ pub fn realize_gram( singular. right now, this causes the Cholesky decomposition to return `None`, leading to a panic when we unrap */ - let base_step_stacked = hess.cholesky().unwrap().solve(&neg_grad_stacked); + let base_step_stacked = hess.clone().cholesky().unwrap().solve(&neg_grad_stacked); let base_step = base_step_stacked.reshape_generic(Dyn(element_dim), Dyn(assembly_dim)); history.base_step.push(base_step.clone()); @@ -269,10 +311,16 @@ pub fn realize_gram( state = better_state; history.backoff_steps.push(backoff_steps); }, - None => return (state.config, false, history) + None => return (state.config, ConfigSubspace::zero(), false, history) }; } - (state.config, state.loss < tol, history) + let success = state.loss < tol; + let tangent = if success { + ConfigSubspace::symmetric_kernel(hess, assembly_dim) + } else { + ConfigSubspace::zero() + }; + (state.config, tangent, success, history) } // --- tests --- @@ -291,7 +339,7 @@ pub mod irisawa { use super::*; - pub fn realize_irisawa_hexlet(scaled_tol: f64) -> (DMatrix, bool, DescentHistory) { + pub fn realize_irisawa_hexlet(scaled_tol: f64) -> (DMatrix, ConfigSubspace, bool, DescentHistory) { let gram = { let mut gram_to_be = PartialMatrix::new(); for s in 0..9 { @@ -399,7 +447,7 @@ mod tests { fn irisawa_hexlet_test() { // solve Irisawa's problem const SCALED_TOL: f64 = 1.0e-12; - let (config, _, _) = realize_irisawa_hexlet(SCALED_TOL); + let (config, _, _, _) = realize_irisawa_hexlet(SCALED_TOL); // check against Irisawa's solution let entry_tol = SCALED_TOL.sqrt(); @@ -409,6 +457,61 @@ mod tests { } } + #[test] + fn tangent_test() { + const SCALED_TOL: f64 = 1.0e-12; + const ELEMENT_DIM: usize = 5; + const ASSEMBLY_DIM: usize = 3; + let gram = { + let mut gram_to_be = PartialMatrix::new(); + for j in 0..3 { + for k in j..3 { + gram_to_be.push_sym(j, k, if j == k { 1.0 } else { -1.0 }); + } + } + gram_to_be + }; + let guess = DMatrix::from_columns(&[ + sphere(0.0, 0.0, 0.0, -2.0), + sphere(0.0, 0.0, 1.0, 1.0), + sphere(0.0, 0.0, -1.0, 1.0) + ]); + let frozen: [_; 5] = std::array::from_fn(|k| (k, 0)); + let (config, tangent, success, history) = realize_gram( + &gram, guess.clone(), &frozen, + SCALED_TOL, 0.5, 0.9, 1.1, 200, 110 + ); + assert_eq!(config, guess); + assert_eq!(success, true); + assert_eq!(history.scaled_loss.len(), 1); + + // confirm that the tangent space has dimension five or less + let ConfigSubspace(ref tangent_basis) = tangent; + assert_eq!(tangent_basis.len(), 5); + + // confirm that the tangent space contains all the motions we expect it + // to. since we've already bounded the dimension of the tangent space, + // this confirms that the tangent space is what we expect it to be + let tangent_motions = vec![ + basis_matrix((0, 1), ELEMENT_DIM, ASSEMBLY_DIM), + basis_matrix((1, 1), ELEMENT_DIM, ASSEMBLY_DIM), + basis_matrix((0, 2), ELEMENT_DIM, ASSEMBLY_DIM), + basis_matrix((1, 2), ELEMENT_DIM, ASSEMBLY_DIM), + DMatrix::::from_column_slice(ELEMENT_DIM, 3, &[ + 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, -1.0, -0.25, -1.0, + 0.0, 0.0, -1.0, 0.25, 1.0 + ]) + ]; + let tol_sq = ((ELEMENT_DIM * ASSEMBLY_DIM) as f64) * SCALED_TOL * SCALED_TOL; + for motion in tangent_motions { + let motion_proj: DMatrix<_> = motion.column_iter().enumerate().map( + |(k, v)| tangent.proj(&v, k) + ).sum(); + assert!((motion - motion_proj).norm_squared() < tol_sq); + } + } + // at the frozen indices, the optimization steps should have exact zeros, // and the realized configuration should match the initial guess #[test] @@ -428,7 +531,7 @@ mod tests { ]); let frozen = [(3, 0), (3, 1)]; println!(); - let (config, success, history) = realize_gram( + let (config, _, success, history) = realize_gram( &gram, guess.clone(), &frozen, 1.0e-12, 0.5, 0.9, 1.1, 200, 110 ); -- 2.34.1 From 7aa69bdfcd58778cdeb8729feb54ecf2eeb646ca Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Sun, 8 Dec 2024 19:59:25 -0800 Subject: [PATCH 2/8] Set the console error panic hook Turn on the browser console panic message output provided by the `console_error_panic_hook` feature. This feature was already enabled by default in our Cargo configuration, but it wasn't actually being used. --- app-proto/src/main.rs | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/app-proto/src/main.rs b/app-proto/src/main.rs index 8a012d3..f961504 100644 --- a/app-proto/src/main.rs +++ b/app-proto/src/main.rs @@ -46,6 +46,10 @@ impl AppState { } fn main() { + // set the console error panic hook + #[cfg(feature = "console_error_panic_hook")] + console_error_panic_hook::set_once(); + sycamore::render(|| { provide_context(AppState::new()); -- 2.34.1 From 58e75871313de7cfc39ccad8d4aba62387a62c8c Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Sun, 8 Dec 2024 20:10:55 -0800 Subject: [PATCH 3/8] Deform the assembly This seems like a good starting point, even though the code is messy and the deformation routine has some numerical quirks. Note that the translation direction is mixed up: the keys are for x, but the velocity field is for z. --- app-proto/src/assembly.rs | 47 ++++++++++++++++++++++++++++++- app-proto/src/display.rs | 58 ++++++++++++++++++++++++++++++++++++--- app-proto/src/engine.rs | 13 +++++++-- 3 files changed, 111 insertions(+), 7 deletions(-) diff --git a/app-proto/src/assembly.rs b/app-proto/src/assembly.rs index 0d63e45..77fbe94 100644 --- a/app-proto/src/assembly.rs +++ b/app-proto/src/assembly.rs @@ -1,4 +1,4 @@ -use nalgebra::{DMatrix, DVector, Vector3}; +use nalgebra::{DMatrix, DVector, DVectorView, Vector3}; use rustc_hash::FxHashMap; use slab::Slab; use std::{collections::BTreeSet, sync::atomic::{AtomicU64, Ordering}}; @@ -152,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 { @@ -266,6 +273,7 @@ impl Assembly { )); 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 @@ -279,4 +287,41 @@ impl Assembly { self.tangent.set_silent(tangent); } } + + // --- deformation --- + + pub fn deform(&self, element_motions: Vec<(ElementKey, DVectorView)>) { + /* 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) { + 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); + + // 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)); + } + + // step the configuration linearly along the tangent space of the + // solution variety + for (_, elt) in self.elements.get_clone_untracked() { + elt.representation.update_silent(|rep| { + let rep_next = &*rep + motion_proj.column(elt.column_index); + rep.set_column(0, &rep_next); + }); + } + + // 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 c39e575..7e1d3a7 100644 --- a/app-proto/src/display.rs +++ b/app-proto/src/display.rs @@ -1,5 +1,5 @@ use core::array; -use nalgebra::{DMatrix, Rotation3, Vector3}; +use nalgebra::{DMatrix, DVector, Rotation3, Vector3}; use sycamore::{prelude::*, motion::create_raf}; use web_sys::{ console, @@ -123,6 +123,10 @@ pub fn Display() -> View { let zoom_out = create_signal(0.0); let turntable = create_signal(false); /* BENCHMARKING */ + // manipulation + let translate_neg_x = create_signal(0.0); + let translate_pos_x = create_signal(0.0); + // change listener let scene_changed = create_signal(true); create_effect(move || { @@ -141,6 +145,7 @@ pub fn Display() -> View { let mut frames_since_last_sample = 0; let mean_frame_interval = create_signal(0.0); + let assembly_for_raf = state.assembly.clone(); on_mount(move || { // timing let mut last_time = 0.0; @@ -153,6 +158,9 @@ pub fn Display() -> View { let mut rotation = DMatrix::::identity(5, 5); let mut location_z: f64 = 5.0; + // manipulation + const TRANSLATION_SPEED: f64 = 0.15; // in length units per second + // display parameters const OPACITY: f32 = 0.5; /* SCAFFOLDING */ const HIGHLIGHT: f32 = 0.2; /* SCAFFOLDING */ @@ -273,6 +281,10 @@ pub fn Display() -> View { let zoom_out_val = zoom_out.get(); let turntable_val = turntable.get(); /* BENCHMARKING */ + // get the manipulation state + let translate_neg_x_val = translate_neg_x.get(); + let translate_pos_x_val = translate_pos_x.get(); + // update the assembly's orientation let ang_vel = { let pitch = pitch_up_val - pitch_down_val; @@ -298,6 +310,30 @@ pub fn Display() -> View { let zoom = zoom_out_val - zoom_in_val; location_z *= (time_step * ZOOM_SPEED * zoom).exp(); + // manipulate the assembly + if state.selection.with(|sel| sel.len() == 1) { + let sel = state.selection.with( + |sel| *sel.into_iter().next().unwrap() + ); + let rep = state.assembly.elements.with_untracked( + |elts| elts[sel].representation.get_clone_untracked() + ); + let vel_field_z = DMatrix::from_column_slice(5, 5, &[ + 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 1.0, + 0.0, 0.0, 2.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0 + ]); + let translate_x = translate_pos_x_val - translate_neg_x_val; + if translate_x != 0.0 { + let vel = translate_x * vel_field_z * rep; + let elt_motion: DVector = time_step * TRANSLATION_SPEED * vel; + assembly_for_raf.deform(vec![(sel, elt_motion.as_view())]); + scene_changed.set(true); + } + } + if scene_changed.get() { /* INSTRUMENTS */ // measure mean frame interval @@ -416,7 +452,7 @@ pub fn Display() -> View { start_animation_loop(); }); - let set_nav_signal = move |event: KeyboardEvent, value: f64| { + let set_nav_signal = move |event: &KeyboardEvent, value: f64| { let mut navigating = true; let shift = event.shift_key(); match event.key().as_str() { @@ -436,6 +472,18 @@ pub fn Display() -> View { } }; + let set_manip_signal = move |event: &KeyboardEvent, value: f64| { + let mut manipulating = true; + match event.key().as_str() { + "d" => translate_pos_x.set(value), + "a" => translate_neg_x.set(value), + _ => manipulating = false + }; + if manipulating { + event.prevent_default(); + } + }; + view! { /* TO DO */ // switch back to integer-valued parameters when that becomes possible @@ -460,7 +508,8 @@ pub fn Display() -> View { turntable.set_fn(|turn| !turn); scene_changed.set(true); } - set_nav_signal(event, 1.0); + set_nav_signal(&event, 1.0); + set_manip_signal(&event, 1.0); } }, on:keyup=move |event: KeyboardEvent| { @@ -474,7 +523,8 @@ pub fn Display() -> View { zoom_in.set(0.0); zoom_out.set(0.0); } else { - set_nav_signal(event, 0.0); + set_nav_signal(&event, 0.0); + set_manip_signal(&event, 0.0); } }, on:blur=move |_| { diff --git a/app-proto/src/engine.rs b/app-proto/src/engine.rs index 6a62579..a48a1af 100644 --- a/app-proto/src/engine.rs +++ b/app-proto/src/engine.rs @@ -87,6 +87,7 @@ impl PartialMatrix { // --- configuration subspaces --- +#[derive(Clone)] pub struct ConfigSubspace(Vec>); impl ConfigSubspace { @@ -99,7 +100,7 @@ impl ConfigSubspace { // of the kernel if its eigenvalue is smaller than the constant `THRESHOLD` fn symmetric_kernel(a: DMatrix, assembly_dim: usize) -> ConfigSubspace { const ELEMENT_DIM: usize = 5; - const THRESHOLD: f64 = 1.0e-9; + const THRESHOLD: f64 = 1.0e-4; let eig = SymmetricEigen::new(a); let eig_vecs = eig.eigenvectors.column_iter(); let eig_pairs = eig.eigenvalues.iter().zip(eig_vecs); @@ -110,15 +111,23 @@ impl ConfigSubspace { ) ) ); + console::log_1(&JsValue::from( + format!("Hessian eigenvalues: {}", eig.eigenvalues) + )); /* DEBUG */ ConfigSubspace(basis.collect()) } + pub fn dim(&self) -> usize { + 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 - fn proj(&self, v: &DVectorView, column_index: usize) -> DMatrix { + 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 -- 2.34.1 From 2906571f326d33f673af3ac4852a91c8b5c67bbe Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Mon, 9 Dec 2024 01:22:54 -0800 Subject: [PATCH 4/8] Correct the translation direction Make the x translation keys translate along the x axis, as intended. --- app-proto/src/display.rs | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/app-proto/src/display.rs b/app-proto/src/display.rs index 7e1d3a7..184b29e 100644 --- a/app-proto/src/display.rs +++ b/app-proto/src/display.rs @@ -318,16 +318,16 @@ pub fn Display() -> View { let rep = state.assembly.elements.with_untracked( |elts| elts[sel].representation.get_clone_untracked() ); - let vel_field_z = DMatrix::from_column_slice(5, 5, &[ - 0.0, 0.0, 0.0, 0.0, 0.0, - 0.0, 0.0, 0.0, 0.0, 0.0, + let vel_field_x = DMatrix::from_column_slice(5, 5, &[ 0.0, 0.0, 0.0, 0.0, 1.0, - 0.0, 0.0, 2.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, + 2.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 ]); let translate_x = translate_pos_x_val - translate_neg_x_val; if translate_x != 0.0 { - let vel = translate_x * vel_field_z * rep; + let vel = translate_x * vel_field_x * rep; let elt_motion: DVector = time_step * TRANSLATION_SPEED * vel; assembly_for_raf.deform(vec![(sel, elt_motion.as_view())]); scene_changed.set(true); -- 2.34.1 From 9f85ce560870fb5a3caff6af1b203da824d88240 Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Mon, 9 Dec 2024 15:58:45 -0800 Subject: [PATCH 5/8] Step elements geodesically instead of linearly This helps prevent small spheres from shrinking during deformations. --- app-proto/src/assembly.rs | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/app-proto/src/assembly.rs b/app-proto/src/assembly.rs index 77fbe94..24f9b93 100644 --- a/app-proto/src/assembly.rs +++ b/app-proto/src/assembly.rs @@ -5,7 +5,7 @@ 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}; +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; @@ -312,12 +312,16 @@ impl Assembly { motion_proj += self.tangent.with(|tan| tan.proj(&elt_motion, column_index)); } - // step the configuration linearly along the tangent space of the - // solution variety + // step each element along the mass shell geodesic that matches its + // velocity in the deformation found above + /* KLUDGE */ + // since our test assemblies only involve 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| { let rep_next = &*rep + motion_proj.column(elt.column_index); - rep.set_column(0, &rep_next); + let normalizer = rep_next.dot(&(&*Q * &rep_next)); + rep.set_column(0, &(rep_next / normalizer)); }); } -- 2.34.1 From 64da1ba5779314eee558b0c6f3fb361ac52e47e2 Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Mon, 9 Dec 2024 21:09:21 -0800 Subject: [PATCH 6/8] Enable translation along all axes --- app-proto/src/display.rs | 55 +++++++++++++++++++++++++++++++--------- 1 file changed, 43 insertions(+), 12 deletions(-) diff --git a/app-proto/src/display.rs b/app-proto/src/display.rs index 184b29e..673985c 100644 --- a/app-proto/src/display.rs +++ b/app-proto/src/display.rs @@ -126,6 +126,10 @@ pub fn Display() -> View { // manipulation let translate_neg_x = create_signal(0.0); let translate_pos_x = create_signal(0.0); + let translate_neg_y = create_signal(0.0); + let translate_pos_y = create_signal(0.0); + let translate_neg_z = create_signal(0.0); + let translate_pos_z = create_signal(0.0); // change listener let scene_changed = create_signal(true); @@ -284,6 +288,10 @@ pub fn Display() -> View { // get the manipulation state let translate_neg_x_val = translate_neg_x.get(); let translate_pos_x_val = translate_pos_x.get(); + let translate_neg_y_val = translate_neg_y.get(); + let translate_pos_y_val = translate_pos_y.get(); + let translate_neg_z_val = translate_neg_z.get(); + let translate_pos_z_val = translate_pos_z.get(); // update the assembly's orientation let ang_vel = { @@ -318,17 +326,21 @@ pub fn Display() -> View { let rep = state.assembly.elements.with_untracked( |elts| elts[sel].representation.get_clone_untracked() ); - let vel_field_x = DMatrix::from_column_slice(5, 5, &[ - 0.0, 0.0, 0.0, 0.0, 1.0, - 0.0, 0.0, 0.0, 0.0, 0.0, - 0.0, 0.0, 0.0, 0.0, 0.0, - 2.0, 0.0, 0.0, 0.0, 0.0, - 0.0, 0.0, 0.0, 0.0, 0.0 - ]); let translate_x = translate_pos_x_val - translate_neg_x_val; - if translate_x != 0.0 { - let vel = translate_x * vel_field_x * rep; - let elt_motion: DVector = time_step * TRANSLATION_SPEED * vel; + let translate_y = translate_pos_y_val - translate_neg_y_val; + let translate_z = translate_pos_z_val - translate_neg_z_val; + if translate_x != 0.0 || translate_y != 0.0 || translate_z != 0.0 { + let vel_field = { + let u = Vector3::new(translate_x, translate_y, translate_z).normalize(); + DMatrix::from_column_slice(5, 5, &[ + 0.0, 0.0, 0.0, 0.0, u[0], + 0.0, 0.0, 0.0, 0.0, u[1], + 0.0, 0.0, 0.0, 0.0, u[2], + 2.0*u[0], 2.0*u[1], 2.0*u[2], 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0 + ]) + }; + let elt_motion: DVector = time_step * TRANSLATION_SPEED * vel_field * rep; assembly_for_raf.deform(vec![(sel, elt_motion.as_view())]); scene_changed.set(true); } @@ -474,9 +486,14 @@ pub fn Display() -> View { let set_manip_signal = move |event: &KeyboardEvent, value: f64| { let mut manipulating = true; + let shift = event.shift_key(); match event.key().as_str() { - "d" => translate_pos_x.set(value), - "a" => translate_neg_x.set(value), + "d" | "D" => translate_pos_x.set(value), + "a" | "A" => translate_neg_x.set(value), + "w" | "W" if shift => translate_neg_z.set(value), + "s" | "S" if shift => translate_pos_z.set(value), + "w" | "W" => translate_pos_y.set(value), + "s" | "S" => translate_neg_y.set(value), _ => manipulating = false }; if manipulating { @@ -495,6 +512,7 @@ pub fn Display() -> View { tabindex="0", on:keydown=move |event: KeyboardEvent| { if event.key() == "Shift" { + // swap navigation inputs roll_cw.set(yaw_right.get()); roll_ccw.set(yaw_left.get()); zoom_in.set(pitch_up.get()); @@ -503,6 +521,12 @@ pub fn Display() -> View { yaw_left.set(0.0); pitch_up.set(0.0); pitch_down.set(0.0); + + // swap manipulation inputs + translate_pos_z.set(translate_neg_y.get()); + translate_neg_z.set(translate_pos_y.get()); + translate_pos_y.set(0.0); + translate_neg_y.set(0.0); } else { if event.key() == "Enter" { /* BENCHMARKING */ turntable.set_fn(|turn| !turn); @@ -514,6 +538,7 @@ pub fn Display() -> View { }, on:keyup=move |event: KeyboardEvent| { if event.key() == "Shift" { + // swap navigation inputs yaw_right.set(roll_cw.get()); yaw_left.set(roll_ccw.get()); pitch_up.set(zoom_in.get()); @@ -522,6 +547,12 @@ pub fn Display() -> View { roll_ccw.set(0.0); zoom_in.set(0.0); zoom_out.set(0.0); + + // swap manipulation inputs + translate_pos_y.set(translate_neg_z.get()); + translate_neg_y.set(translate_pos_z.get()); + translate_pos_z.set(0.0); + translate_neg_z.set(0.0); } else { set_nav_signal(&event, 0.0); set_manip_signal(&event, 0.0); -- 2.34.1 From c87367a276d7bbb6cb7280d058da013c1e060888 Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Tue, 10 Dec 2024 01:56:10 -0800 Subject: [PATCH 7/8] Tweak comment wording --- app-proto/src/assembly.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/app-proto/src/assembly.rs b/app-proto/src/assembly.rs index 24f9b93..1684716 100644 --- a/app-proto/src/assembly.rs +++ b/app-proto/src/assembly.rs @@ -315,7 +315,7 @@ impl Assembly { // step each element along the mass shell geodesic that matches its // velocity in the deformation found above /* KLUDGE */ - // since our test assemblies only involve spheres, we assume that every + // 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| { -- 2.34.1 From 90834fbb937fdce1acf7861fdb64d36c0d2021f4 Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Wed, 11 Dec 2024 12:59:57 -0800 Subject: [PATCH 8/8] Adapt `symmetric_kernel` for non-WASM targets The examples call `engine::realize_gram`, which now includes a call to `symmetric_kernel`, so we need to make sure that `symmetric_kernel` can run on whatever target Cargo uses for examples. For that target on my machine, `console::log_1` panics with the message "function not implemented on non-`wasm32` targets". --- app-proto/src/engine.rs | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/app-proto/src/engine.rs b/app-proto/src/engine.rs index a48a1af..2f06035 100644 --- a/app-proto/src/engine.rs +++ b/app-proto/src/engine.rs @@ -111,9 +111,14 @@ impl ConfigSubspace { ) ) ); + + /* DEBUG */ + // print the eigenvalues + #[cfg(all(target_family = "wasm", target_os = "unknown"))] console::log_1(&JsValue::from( - format!("Hessian eigenvalues: {}", eig.eigenvalues) - )); /* DEBUG */ + format!("Eigenvalues used to find kernel: {}", eig.eigenvalues) + )); + ConfigSubspace(basis.collect()) } -- 2.34.1