From 51df7defc52a0efc93f6a805187859aa28c1e51d Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Thu, 23 Jan 2025 12:16:07 -0800 Subject: [PATCH] Test the kaleidocycle's first-order motions --- app-proto/src/engine.rs | 144 +++++++++++++++++++++++++++++++++++++--- 1 file changed, 136 insertions(+), 8 deletions(-) diff --git a/app-proto/src/engine.rs b/app-proto/src/engine.rs index 0756fac..57243b8 100644 --- a/app-proto/src/engine.rs +++ b/app-proto/src/engine.rs @@ -491,7 +491,7 @@ pub mod irisawa { #[cfg(test)] mod tests { use nalgebra::Vector3; - use std::f64::consts::FRAC_1_SQRT_2; + use std::{array, f64::consts::{FRAC_1_SQRT_2, PI}, iter}; use super::{*, irisawa::realize_irisawa_hexlet}; @@ -555,7 +555,7 @@ mod tests { } #[test] - fn tangent_test() { + fn tangent_test_three_spheres() { const SCALED_TOL: f64 = 1.0e-12; let gram = { let mut gram_to_be = PartialMatrix::new(); @@ -620,14 +620,142 @@ mod tests { } } - fn translation(u: Vector3) -> DMatrix { + fn translation_motion_unif(vel: &Vector3, assembly_dim: usize) -> Vec> { + let mut elt_motion = DVector::zeros(4); + elt_motion.fixed_rows_mut::<3>(0).copy_from(vel); + iter::repeat(elt_motion).take(assembly_dim).collect() + } + + fn rotation_motion_unif(ang_vel: &Vector3, points: Vec>) -> Vec> { + points.into_iter().map( + |pt| { + let vel = ang_vel.cross(&pt.fixed_rows::<3>(0)); + let mut elt_motion = DVector::zeros(4); + elt_motion.fixed_rows_mut::<3>(0).copy_from(&vel); + elt_motion + } + ).collect() + } + + #[test] + fn tangent_test_kaleidocycle() { + const SCALED_TOL: f64 = 1.0e-12; + // set up a kaleidocycle, made of points with fixed distances between + // them, and find its tangent space + const N_POINTS: usize = 12; + const N_HINGES: usize = 6; + 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 = { + 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.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); + + // list some motions that should form a basis for the tangent space of + // the solution variety + let element_dim = guess.nrows(); + let assembly_dim = guess.ncols(); + let tangent_motions_unif = vec![ + // the translations along the coordinate axes + translation_motion_unif(&Vector3::new(1.0, 0.0, 0.0), assembly_dim), + translation_motion_unif(&Vector3::new(0.0, 1.0, 0.0), assembly_dim), + translation_motion_unif(&Vector3::new(0.0, 0.0, 1.0), assembly_dim), + + // the rotations about the coordinate axes + rotation_motion_unif(&Vector3::new(1.0, 0.0, 0.0), guess.column_iter().collect()), + rotation_motion_unif(&Vector3::new(0.0, 1.0, 0.0), guess.column_iter().collect()), + rotation_motion_unif(&Vector3::new(0.0, 0.0, 1.0), guess.column_iter().collect()), + + // the twist motion. more precisely: a motion that keeps the center + // of mass stationary and preserves the distances between the + // vertices to first order. this has to be the twist as long as: + // - twisting is the kaleidocycle's only internal degree of + // freedom + // - every first-order motion of the kaleidocycle comes from an + // actual motion + (0..N_HINGES).step_by(2).flat_map( + |n| { + let ang_vert = ((n + 1) as f64) * PI/3.0; + let vel_vert_x = 4.0 * ang_vert.cos(); + let vel_vert_y = 4.0 * ang_vert.sin(); + [ + DVector::from_column_slice(&[0.0, 0.0, 5.0, 0.0]), + DVector::from_column_slice(&[0.0, 0.0, 1.0, 0.0]), + DVector::from_column_slice(&[-vel_vert_x, -vel_vert_y, -3.0, 0.0]), + DVector::from_column_slice(&[vel_vert_x, vel_vert_y, -3.0, 0.0]) + ] + } + ).collect::>() + ]; + let tangent_motions_std = tangent_motions_unif.iter().map( + |motion| DMatrix::from_columns( + &guess.column_iter().zip(motion).map( + |(v, elt_motion)| local_unif_to_std(v) * elt_motion + ).collect::>() + ) + ).collect::>(); + + // confirm that the dimension of the tangent space is no greater than + // expected + assert_eq!(tangent.basis_std.len(), tangent_motions_unif.len()); + + // 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 tol_sq = ((element_dim * assembly_dim) as f64) * SCALED_TOL * SCALED_TOL; + for (motion_unif, motion_std) in tangent_motions_unif.into_iter().zip(tangent_motions_std) { + let motion_proj: DMatrix<_> = motion_unif.into_iter().enumerate().map( + |(k, v)| tangent.proj(&v.as_view(), k) + ).sum(); + assert!((motion_std - motion_proj).norm_squared() < tol_sq); + } + } + + fn translation(dis: Vector3) -> DMatrix { const ELEMENT_DIM: usize = 5; DMatrix::from_column_slice(ELEMENT_DIM, ELEMENT_DIM, &[ - 1.0, 0.0, 0.0, 0.0, u[0], - 0.0, 1.0, 0.0, 0.0, u[1], - 0.0, 0.0, 1.0, 0.0, u[2], - 2.0*u[0], 2.0*u[1], 2.0*u[2], 1.0, u.norm_squared(), - 0.0, 0.0, 0.0, 0.0, 1.0 + 1.0, 0.0, 0.0, 0.0, dis[0], + 0.0, 1.0, 0.0, 0.0, dis[1], + 0.0, 0.0, 1.0, 0.0, dis[2], + 2.0*dis[0], 2.0*dis[1], 2.0*dis[2], 1.0, dis.norm_squared(), + 0.0, 0.0, 0.0, 0.0, 1.0 ]) }