From 03da831c9a73375c3bf9e778b806d40e7aa70b89 Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Fri, 17 Jan 2025 00:58:35 -0800 Subject: [PATCH 1/4] Try projecting with the uniform inner product This projection should be invariant under Euclidean motions. For now, we assume that all of the elements are spheres. --- app-proto/src/assembly.rs | 12 +++++- app-proto/src/display.rs | 38 ++++++++++++------- app-proto/src/engine.rs | 80 ++++++++++++++++++++++++++++++--------- 3 files changed, 97 insertions(+), 33 deletions(-) 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) }; -- 2.34.1 From a05a2e1d5498f1d2263d3622133cece4d748773d Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Mon, 20 Jan 2025 15:44:20 -0800 Subject: [PATCH 2/4] Implement the uniform inner product for points --- app-proto/examples/kaleidocycle.rs | 81 ++++++++++++++++++++++++++++++ app-proto/run-examples | 1 + app-proto/src/engine.rs | 40 +++++++++++---- 3 files changed, 112 insertions(+), 10 deletions(-) create mode 100644 app-proto/examples/kaleidocycle.rs diff --git a/app-proto/examples/kaleidocycle.rs b/app-proto/examples/kaleidocycle.rs new file mode 100644 index 0000000..effb226 --- /dev/null +++ b/app-proto/examples/kaleidocycle.rs @@ -0,0 +1,81 @@ +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 bc7e933..52173b0 100755 --- a/app-proto/run-examples +++ b/app-proto/run-examples @@ -9,3 +9,4 @@ 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/engine.rs b/app-proto/src/engine.rs index c6b6150..9db6db0 100644 --- a/app-proto/src/engine.rs +++ b/app-proto/src/engine.rs @@ -132,6 +132,9 @@ impl ConfigSubspace { // 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) )); @@ -236,19 +239,36 @@ 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] - ]) + 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 -- 2.34.1 From 01261aed91d4b4dae8f898b904a29720ddb5cf9a Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Tue, 21 Jan 2025 18:04:21 -0800 Subject: [PATCH 3/4] Write kaleidocycle nudge in uniform coordinates In the previous commit, the nudge was written in standard coordinates, so the example shouldn't have worked. However, by sheer dumb luck, the standard and uniform coordinates match for this particular nudge. In fact, the expression used before to get the standard coordinates may even produce equivalent floating point values as the expression used here to get the uniform coordinates. That would explain why the example prints exactly the same output in this commit. --- app-proto/examples/kaleidocycle.rs | 23 +++++++---------------- 1 file changed, 7 insertions(+), 16 deletions(-) diff --git a/app-proto/examples/kaleidocycle.rs b/app-proto/examples/kaleidocycle.rs index effb226..3448b87 100644 --- a/app-proto/examples/kaleidocycle.rs +++ b/app-proto/examples/kaleidocycle.rs @@ -1,4 +1,4 @@ -use nalgebra::DMatrix; +use nalgebra::{DMatrix, DVector}; use std::{array, f64::consts::PI}; use dyna3::engine::{Q, point, realize_gram, PartialMatrix}; @@ -59,22 +59,13 @@ 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, 0.0]); + let down = -&up; 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) - ] - } + |n| [ + tangent.proj(&up.as_view(), n), + tangent.proj(&down.as_view(), n+1) + ] ).sum(); let normalization = 5.0 / twist_motion[(2, 0)]; print!("Twist motion:{}", normalization * twist_motion); -- 2.34.1 From 21cefa9f8a2f297b6cf36fbc31ae5249cca3a88f Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Tue, 21 Jan 2025 18:52:01 -0800 Subject: [PATCH 4/4] Do `symmetric_kernel` in projection coordinates Instead of finding the kernel in the standard coordinates and then expressing it in the projection coordinates, work in the projection coordinates from the beginning by applying a change of basis to the input matrix. --- app-proto/examples/kaleidocycle.rs | 2 +- app-proto/src/display.rs | 2 +- app-proto/src/engine.rs | 46 +++++++++++++++--------------- 3 files changed, 25 insertions(+), 25 deletions(-) diff --git a/app-proto/examples/kaleidocycle.rs b/app-proto/examples/kaleidocycle.rs index 3448b87..88116d3 100644 --- a/app-proto/examples/kaleidocycle.rs +++ b/app-proto/examples/kaleidocycle.rs @@ -59,7 +59,7 @@ 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, 0.0]); + 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| [ diff --git a/app-proto/src/display.rs b/app-proto/src/display.rs index bbfe059..4e0c7e4 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, 0.0] + &[u[0], u[1], u[2], SHRINKING_SPEED * shrink] ) }; assembly_for_raf.deform( diff --git a/app-proto/src/engine.rs b/app-proto/src/engine.rs index 9db6db0..0f1fd36 100644 --- a/app-proto/src/engine.rs +++ b/app-proto/src/engine.rs @@ -107,13 +107,14 @@ 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, expressed in the standard coordinates - const ELEMENT_DIM: usize = 5; - const THRESHOLD: f64 = 1.0e-4; - let eig = SymmetricEigen::new(a); + // 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); let eig_vecs = eig.eigenvectors.column_iter(); let eig_pairs = eig.eigenvalues.iter().zip(eig_vecs); - let basis_std = DMatrix::from_columns( + let basis_proj = DMatrix::from_columns( eig_pairs.filter_map( |(λ, v)| (λ.abs() < THRESHOLD).then_some(v) ).collect::>().as_slice() @@ -126,29 +127,27 @@ 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; + // express the basis in the standard coordinates + let basis_std = proj_to_std * &basis_proj; // 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) + format!("Basis in projection coordinates:{}", basis_proj) )); + const ELEMENT_DIM: usize = 5; + const UNIFORM_DIM: usize = 4; ConfigSubspace { assembly_dim: assembly_dim, - basis_std: basis_std_orth.column_iter().map( + basis_std: basis_std.column_iter().map( |v| Into::>::into( v.reshape_generic(Dyn(ELEMENT_DIM), Dyn(assembly_dim)) ) ).collect(), - basis_proj: basis_proj_orth.column_iter().map( + basis_proj: basis_proj.column_iter().map( |v| Into::>::into( - v.reshape_generic(Dyn(ELEMENT_DIM), Dyn(assembly_dim)) + v.reshape_generic(Dyn(UNIFORM_DIM), Dyn(assembly_dim)) ) ).collect() } @@ -247,26 +246,25 @@ 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, ELEMENT_DIM, &[ + DMatrix::from_column_slice(ELEMENT_DIM, UNIFORM_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, &[ + DMatrix::from_column_slice(ELEMENT_DIM, UNIFORM_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] + curv*v[0], curv*v[1], curv*v[2], curv*v[3], curv*v[4] + 1.0 ]) } } @@ -401,11 +399,13 @@ 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); + 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); for n in 0..assembly_dim { - let block_start = element_dim * n; + let block_start = (element_dim * n, UNIFORM_DIM * n); unif_to_std - .view_mut((block_start, block_start), (element_dim, element_dim)) + .view_mut(block_start, (element_dim, UNIFORM_DIM)) .copy_from(&local_unif_to_std(state.config.column(n))); } -- 2.34.1