From 58e75871313de7cfc39ccad8d4aba62387a62c8c Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Sun, 8 Dec 2024 20:10:55 -0800 Subject: [PATCH] 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