diff --git a/app-proto/examples/kaleidocycle.rs b/app-proto/examples/kaleidocycle.rs index 88116d3..effb226 100644 --- a/app-proto/examples/kaleidocycle.rs +++ b/app-proto/examples/kaleidocycle.rs @@ -1,4 +1,4 @@ -use nalgebra::{DMatrix, DVector}; +use nalgebra::DMatrix; use std::{array, f64::consts::PI}; use dyna3::engine::{Q, point, realize_gram, PartialMatrix}; @@ -59,13 +59,22 @@ fn main() { println!("Loss: {}\n", history.scaled_loss.last().unwrap()); // find the kaleidocycle's twist motion - let up = DVector::from_column_slice(&[0.0, 0.0, 1.0, 0.0]); - let down = -&up; let twist_motion: DMatrix<_> = (0..N_POINTS).step_by(4).flat_map( - |n| [ - tangent.proj(&up.as_view(), n), - tangent.proj(&down.as_view(), n+1) - ] + |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); diff --git a/app-proto/src/display.rs b/app-proto/src/display.rs index 4e0c7e4..bbfe059 100644 --- a/app-proto/src/display.rs +++ b/app-proto/src/display.rs @@ -346,7 +346,7 @@ pub fn Display() -> View { Vector3::zeros() }; time_step * DVector::from_column_slice( - &[u[0], u[1], u[2], SHRINKING_SPEED * shrink] + &[u[0], u[1], u[2], SHRINKING_SPEED * shrink, 0.0] ) }; assembly_for_raf.deform( diff --git a/app-proto/src/engine.rs b/app-proto/src/engine.rs index 0f1fd36..9db6db0 100644 --- a/app-proto/src/engine.rs +++ b/app-proto/src/engine.rs @@ -107,14 +107,13 @@ impl ConfigSubspace { // 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. the basis is expressed in the projection - // coordinates, and it's orthonormal with respect to the projection - // inner product - const THRESHOLD: f64 = 0.1; - let eig = SymmetricEigen::new(proj_to_std.tr_mul(&a) * &proj_to_std); + // 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_proj = DMatrix::from_columns( + let basis_std = DMatrix::from_columns( eig_pairs.filter_map( |(λ, v)| (λ.abs() < THRESHOLD).then_some(v) ).collect::>().as_slice() @@ -127,27 +126,29 @@ impl ConfigSubspace { format!("Eigenvalues used to find kernel:{}", eig.eigenvalues) )); - // express the basis in the standard coordinates - let basis_std = proj_to_std * &basis_proj; + // 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) + format!("Basis in projection coordinates:{}", basis_proj_orth) )); - const ELEMENT_DIM: usize = 5; - const UNIFORM_DIM: usize = 4; ConfigSubspace { assembly_dim: assembly_dim, - basis_std: basis_std.column_iter().map( + 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.column_iter().map( + basis_proj: basis_proj_orth.column_iter().map( |v| Into::>::into( - v.reshape_generic(Dyn(UNIFORM_DIM), Dyn(assembly_dim)) + v.reshape_generic(Dyn(ELEMENT_DIM), Dyn(assembly_dim)) ) ).collect() } @@ -246,25 +247,26 @@ fn basis_matrix(index: (usize, usize), nrows: usize, ncols: usize) -> DMatrix) -> DMatrix { const ELEMENT_DIM: usize = 5; - const UNIFORM_DIM: usize = 4; 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, UNIFORM_DIM, &[ + 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, UNIFORM_DIM, &[ + 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 + 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] ]) } } @@ -399,13 +401,11 @@ pub fn realize_gram( let success = state.loss < tol; let tangent = if success { // express the uniform basis in the standard basis - const UNIFORM_DIM: usize = 4; - let total_dim_unif = UNIFORM_DIM * assembly_dim; - let mut unif_to_std = DMatrix::::zeros(total_dim, total_dim_unif); + let mut unif_to_std = DMatrix::::zeros(total_dim, total_dim); for n in 0..assembly_dim { - let block_start = (element_dim * n, UNIFORM_DIM * n); + let block_start = element_dim * n; unif_to_std - .view_mut(block_start, (element_dim, UNIFORM_DIM)) + .view_mut((block_start, block_start), (element_dim, element_dim)) .copy_from(&local_unif_to_std(state.config.column(n))); }