glen/dyna3
1
0
Fork 0
Code Issues 24 Pull Requests Actions Projects Releases 2 Wiki Activity

Test the kaleidocycle's first-order motions

Browse Source
This commit is contained in:
Aaron Fenyes 2025-01-23 12:16:07 -08:00
parent 1e3505dd01
commit 51df7defc5
1 changed files with 136 additions and 8 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

142
app-proto/src/engine.rs
View File

@ -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,13 +620,141 @@ mod tests {
}
}
fn translation(u: Vector3<f64>) -> DMatrix<f64> {
fn translation_motion_unif(vel: &Vector3<f64>, assembly_dim: usize) -> Vec<DVector<f64>> {
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<f64>, points: Vec<DVectorView<f64>>) -> Vec<DVector<f64>> {
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::<Vec<_>>();
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::<Vec<_>>()
];
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::<Vec<_>>()
)
).collect::<Vec<_>>();
// 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<f64>) -> DMatrix<f64> {
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(),
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
])
}