diff --git a/app-proto/examples/kaleidocycle.rs b/app-proto/examples/kaleidocycle.rs deleted file mode 100644 index effb226..0000000 --- a/app-proto/examples/kaleidocycle.rs +++ /dev/null @@ -1,81 +0,0 @@ -use nalgebra::DMatrix; -use std::{array, f64::consts::PI}; - -use dyna3::engine::{Q, point, realize_gram, PartialMatrix}; - -fn main() { - // set up a kaleidocycle, made of points with fixed distances between them, - // and find its tangent space - const N_POINTS: usize = 12; - let gram = { - let mut gram_to_be = PartialMatrix::new(); - for block in (0..N_POINTS).step_by(2) { - let block_next = (block + 2) % N_POINTS; - for j in 0..2 { - // diagonal and hinge edges - for k in j..2 { - gram_to_be.push_sym(block + j, block + k, if j == k { 0.0 } else { -0.5 }); - } - - // non-hinge edges - for k in 0..2 { - gram_to_be.push_sym(block + j, block_next + k, -0.625); - } - } - } - gram_to_be - }; - let guess = { - const N_HINGES: usize = 6; - let guess_elts = (0..N_HINGES).step_by(2).flat_map( - |n| { - let ang_hor = (n as f64) * PI/3.0; - let ang_vert = ((n + 1) as f64) * PI/3.0; - let x_vert = ang_vert.cos(); - let y_vert = ang_vert.sin(); - [ - point(0.0, 0.0, 0.0), - point(ang_hor.cos(), ang_hor.sin(), 0.0), - point(x_vert, y_vert, -0.5), - point(x_vert, y_vert, 0.5) - ] - } - ).collect::>(); - DMatrix::from_columns(&guess_elts) - }; - let frozen: [_; N_POINTS] = array::from_fn(|k| (3, k)); - let (config, tangent, success, history) = realize_gram( - &gram, guess, &frozen, - 1.0e-12, 0.5, 0.9, 1.1, 200, 110 - ); - print!("Completed Gram matrix:{}", config.tr_mul(&*Q) * &config); - print!("Configuration:{}", config); - if success { - println!("Target accuracy achieved!"); - } else { - println!("Failed to reach target accuracy"); - } - println!("Steps: {}", history.scaled_loss.len() - 1); - println!("Loss: {}\n", history.scaled_loss.last().unwrap()); - - // find the kaleidocycle's twist motion - let twist_motion: DMatrix<_> = (0..N_POINTS).step_by(4).flat_map( - |n| { - let up_field = { - 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 - ]) - }; - [ - tangent.proj(&(&up_field * config.column(n)).as_view(), n), - tangent.proj(&(-&up_field * config.column(n+1)).as_view(), n+1) - ] - } - ).sum(); - let normalization = 5.0 / twist_motion[(2, 0)]; - print!("Twist motion:{}", normalization * twist_motion); -} \ No newline at end of file diff --git a/app-proto/run-examples b/app-proto/run-examples index 52173b0..bc7e933 100755 --- a/app-proto/run-examples +++ b/app-proto/run-examples @@ -9,4 +9,3 @@ cargo run --example irisawa-hexlet cargo run --example three-spheres cargo run --example point-on-sphere -cargo run --example kaleidocycle \ No newline at end of file diff --git a/app-proto/src/assembly.rs b/app-proto/src/assembly.rs index 9b0e065..278a8f9 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, local_unif_to_std, ConfigSubspace, PartialMatrix, Q}; +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; @@ -120,7 +120,6 @@ 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> @@ -360,14 +359,7 @@ 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); - 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; + target_column += elt_motion.velocity; } } diff --git a/app-proto/src/display.rs b/app-proto/src/display.rs index bbfe059..d8ee23c 100644 --- a/app-proto/src/display.rs +++ b/app-proto/src/display.rs @@ -130,8 +130,6 @@ 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); @@ -166,7 +164,6 @@ 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 */ @@ -295,8 +292,6 @@ 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 = { @@ -328,27 +323,24 @@ 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; - 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] - ) + 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![ ElementMotion { @@ -509,8 +501,6 @@ 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 9db6db0..6a0202e 100644 --- a/app-proto/src/engine.rs +++ b/app-proto/src/engine.rs @@ -90,33 +90,32 @@ impl PartialMatrix { #[derive(Clone)] pub struct ConfigSubspace { assembly_dim: usize, - basis_std: Vec>, - basis_proj: Vec> + basis: Vec> } impl ConfigSubspace { pub fn zero(assembly_dim: usize) -> ConfigSubspace { ConfigSubspace { assembly_dim: assembly_dim, - basis_proj: Vec::new(), - basis_std: Vec::new() + basis: 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, proj_to_std: DMatrix, assembly_dim: usize) -> ConfigSubspace { - // find a basis for the kernel, expressed in the standard coordinates + fn symmetric_kernel(a: DMatrix, assembly_dim: usize) -> ConfigSubspace { 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_std = DMatrix::from_columns( - eig_pairs.filter_map( - |(λ, v)| (λ.abs() < THRESHOLD).then_some(v) - ).collect::>().as_slice() + let basis = eig_pairs.filter_map( + |(λ, v)| (λ.abs() < THRESHOLD).then_some( + Into::>::into( + v.reshape_generic(Dyn(ELEMENT_DIM), Dyn(assembly_dim)) + ) + ) ); /* DEBUG */ @@ -126,53 +125,30 @@ 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; - - // print the projection basis in projection coordinates - #[cfg(all(target_family = "wasm", target_os = "unknown"))] - console::log_1(&JsValue::from( - format!("Basis in projection coordinates:{}", basis_proj_orth) - )); - ConfigSubspace { assembly_dim: assembly_dim, - 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() + basis: basis.collect() } } pub fn dim(&self) -> usize { - self.basis_std.len() + self.basis.len() } pub fn assembly_dim(&self) -> usize { self.assembly_dim } - // 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 + // 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` 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_proj.iter().zip(self.basis_std.iter()).map( - |(b_proj, b_std)| b_proj.column(column_index).dot(&v) * b_std + self.basis.iter().map( + |b| b.column(column_index).dot(&v) * b ).sum() } } @@ -239,38 +215,6 @@ fn basis_matrix(index: (usize, usize), nrows: usize, ncols: usize) -> DMatrix) -> DMatrix { - const ELEMENT_DIM: usize = 5; - let curv = 2.0*v[3]; - if v.dot(&(&*Q * v)) < 0.5 { - // `v` represents a point. the normalization condition says that the - // curvature component of `v` is 1/2 - 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], - v[0], v[1], v[2], v[3], v[4], - 0.0, 0.0, 0.0, 0.0, 1.0 - ]) - } else { - // `v` represents a sphere. the normalization condition says that the - // Lorentz product of `v` with itself is 1 - 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, @@ -400,17 +344,7 @@ pub fn realize_gram( } let success = state.loss < tol; let tangent = if success { - // 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) + ConfigSubspace::symmetric_kernel(hess, assembly_dim) } else { ConfigSubspace::zero(assembly_dim) };