diff --git a/app-proto/src/assembly.rs b/app-proto/src/assembly.rs index 278a8f9..9b0e065 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, Q}; +use crate::engine::{realize_gram, local_unif_to_std, ConfigSubspace, PartialMatrix, Q}; // the types of the keys we use to access an assembly's elements and constraints pub type ElementKey = usize; @@ -120,6 +120,7 @@ pub struct Constraint { pub active: Signal } +// the velocity is expressed in uniform coordinates pub struct ElementMotion<'a> { pub key: ElementKey, pub velocity: DVectorView<'a, f64> @@ -359,7 +360,14 @@ impl Assembly { // 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; + let unif_to_std = self.elements.with_untracked( + |elts| { + elts[elt_motion.key].representation.with_untracked( + |rep| local_unif_to_std(rep.as_view()) + ) + } + ); + target_column += unif_to_std * elt_motion.velocity; } } diff --git a/app-proto/src/display.rs b/app-proto/src/display.rs index d8ee23c..bbfe059 100644 --- a/app-proto/src/display.rs +++ b/app-proto/src/display.rs @@ -130,6 +130,8 @@ pub fn Display() -> View { let translate_pos_y = create_signal(0.0); let translate_neg_z = create_signal(0.0); let translate_pos_z = create_signal(0.0); + let shrink_neg = create_signal(0.0); + let shrink_pos = create_signal(0.0); // change listener let scene_changed = create_signal(true); @@ -164,6 +166,7 @@ pub fn Display() -> View { // manipulation const TRANSLATION_SPEED: f64 = 0.15; // in length units per second + const SHRINKING_SPEED: f64 = 0.15; // in length units per second // display parameters const OPACITY: f32 = 0.5; /* SCAFFOLDING */ @@ -292,6 +295,8 @@ pub fn Display() -> View { 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(); + let shrink_neg_val = shrink_neg.get(); + let shrink_pos_val = shrink_pos.get(); // update the assembly's orientation let ang_vel = { @@ -323,24 +328,27 @@ pub fn Display() -> View { 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 translate_x = translate_pos_x_val - translate_neg_x_val; 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 shrink = shrink_pos_val - shrink_neg_val; + let translating = + translate_x != 0.0 + || translate_y != 0.0 + || translate_z != 0.0; + if translating || shrink != 0.0 { + let elt_motion = { + let u = if translating { + TRANSLATION_SPEED * Vector3::new( + translate_x, translate_y, translate_z + ).normalize() + } else { + Vector3::zeros() + }; + time_step * DVector::from_column_slice( + &[u[0], u[1], u[2], SHRINKING_SPEED * shrink, 0.0] + ) }; - let elt_motion: DVector = time_step * TRANSLATION_SPEED * vel_field * rep; assembly_for_raf.deform( vec![ ElementMotion { @@ -501,6 +509,8 @@ pub fn Display() -> View { "s" | "S" if shift => translate_pos_z.set(value), "w" | "W" => translate_pos_y.set(value), "s" | "S" => translate_neg_y.set(value), + "]" | "}" => shrink_neg.set(value), + "[" | "{" => shrink_pos.set(value), _ => manipulating = false }; if manipulating { diff --git a/app-proto/src/engine.rs b/app-proto/src/engine.rs index 6a0202e..c6b6150 100644 --- a/app-proto/src/engine.rs +++ b/app-proto/src/engine.rs @@ -90,32 +90,33 @@ impl PartialMatrix { #[derive(Clone)] pub struct ConfigSubspace { assembly_dim: usize, - basis: Vec> + basis_std: Vec>, + basis_proj: Vec> } impl ConfigSubspace { pub fn zero(assembly_dim: usize) -> ConfigSubspace { ConfigSubspace { assembly_dim: assembly_dim, - basis: Vec::new() + basis_proj: Vec::new(), + basis_std: 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 { + fn symmetric_kernel(a: DMatrix, proj_to_std: DMatrix, assembly_dim: usize) -> ConfigSubspace { + // find a basis for the kernel, expressed in the standard coordinates const ELEMENT_DIM: usize = 5; 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); - let basis = eig_pairs.filter_map( - |(λ, v)| (λ.abs() < THRESHOLD).then_some( - Into::>::into( - v.reshape_generic(Dyn(ELEMENT_DIM), Dyn(assembly_dim)) - ) - ) + let basis_std = DMatrix::from_columns( + eig_pairs.filter_map( + |(λ, v)| (λ.abs() < THRESHOLD).then_some(v) + ).collect::>().as_slice() ); /* DEBUG */ @@ -125,30 +126,50 @@ impl ConfigSubspace { format!("Eigenvalues used to find kernel:{}", eig.eigenvalues) )); + // express the basis in the projection coordinates + let basis_proj = proj_to_std.clone().qr().solve(&basis_std).unwrap(); + + // orthonormalize the basis with respect to the projection inner product + let basis_proj_orth = basis_proj.qr().q(); + let basis_std_orth = proj_to_std * &basis_proj_orth; + console::log_1(&JsValue::from( + format!("Basis in projection coordinates:{}", basis_proj_orth) + )); + ConfigSubspace { assembly_dim: assembly_dim, - basis: basis.collect() + basis_std: basis_std_orth.column_iter().map( + |v| Into::>::into( + v.reshape_generic(Dyn(ELEMENT_DIM), Dyn(assembly_dim)) + ) + ).collect(), + basis_proj: basis_proj_orth.column_iter().map( + |v| Into::>::into( + v.reshape_generic(Dyn(ELEMENT_DIM), Dyn(assembly_dim)) + ) + ).collect() } } pub fn dim(&self) -> usize { - self.basis.len() + self.basis_std.len() } pub fn assembly_dim(&self) -> usize { self.assembly_dim } - // find the projection onto this subspace, with respect to the Euclidean - // inner product, of the motion where the element with the given column - // index has velocity `v` + // find the projection onto this subspace of the motion where the element + // with the given column index has velocity `v`. the velocity is given in + // projection coordinates, and the projection is done with respect to the + // projection inner product 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 + self.basis_proj.iter().zip(self.basis_std.iter()).map( + |(b_proj, b_std)| b_proj.column(column_index).dot(&v) * b_std ).sum() } } @@ -215,6 +236,21 @@ fn basis_matrix(index: (usize, usize), nrows: usize, ncols: usize) -> DMatrix) -> DMatrix { + const ELEMENT_DIM: usize = 5; + let curv = 2.0*v[3]; + DMatrix::from_column_slice(ELEMENT_DIM, ELEMENT_DIM, &[ + curv, 0.0, 0.0, 0.0, v[0], + 0.0, curv, 0.0, 0.0, v[1], + 0.0, 0.0, curv, 0.0, v[2], + curv*v[0], curv*v[1], curv*v[2], curv*v[3], curv*v[4] + 1.0, + v[0], v[1], v[2], v[3], v[4] + ]) +} + // use backtracking line search to find a better configuration fn seek_better_config( gram: &PartialMatrix, @@ -344,7 +380,17 @@ pub fn realize_gram( } let success = state.loss < tol; let tangent = if success { - ConfigSubspace::symmetric_kernel(hess, assembly_dim) + // express the uniform basis in the standard basis + let mut unif_to_std = DMatrix::::zeros(total_dim, total_dim); + for n in 0..assembly_dim { + let block_start = element_dim * n; + unif_to_std + .view_mut((block_start, block_start), (element_dim, element_dim)) + .copy_from(&local_unif_to_std(state.config.column(n))); + } + + // find the kernel of the Hessian. give it the uniform inner product + ConfigSubspace::symmetric_kernel(hess, unif_to_std, assembly_dim) } else { ConfigSubspace::zero(assembly_dim) };