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

@@ -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<f64>)>) {
+        /* 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();
+    }
 }

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::<f64>::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<f64> = 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 |_| {

app-proto/src/engine.rs

@@ -87,6 +87,7 @@ impl PartialMatrix {
 
 // --- configuration subspaces ---
 
+#[derive(Clone)]
 pub struct ConfigSubspace(Vec<DMatrix<f64>>);
 
 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<f64>, 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<f64>, column_index: usize) -> DMatrix<f64> {
+    pub fn proj(&self, v: &DVectorView<f64>, column_index: usize) -> DMatrix<f64> {
         let ConfigSubspace(basis) = self;
         basis.into_iter().map(
             |b| b.column(column_index).dot(&v) * b