Compare commits

..

No commits in common. "6df0e855cf3a7ca144c2388965215ad6888a8bcb" and "90834fbb937fdce1acf7861fdb64d36c0d2021f4" have entirely different histories.

3 changed files with 52 additions and 134 deletions

View file

@ -33,9 +33,8 @@ pub struct Element {
pub serial: u64, pub serial: u64,
// the configuration matrix column index that was assigned to this element // the configuration matrix column index that was assigned to this element
// last time the assembly was realized, or `None` if the element has never // last time the assembly was realized
// been through a realization column_index: usize
column_index: Option<usize>
} }
impl Element { impl Element {
@ -63,7 +62,7 @@ impl Element {
representation: create_signal(representation), representation: create_signal(representation),
constraints: create_signal(BTreeSet::default()), constraints: create_signal(BTreeSet::default()),
serial: serial, serial: serial,
column_index: None column_index: 0
} }
} }
@ -120,13 +119,6 @@ pub struct Constraint {
pub active: Signal<bool> pub active: Signal<bool>
} }
pub struct ElementMotion<'a> {
pub key: ElementKey,
pub velocity: DVectorView<'a, f64>
}
type AssemblyMotion<'a> = Vec<ElementMotion<'a>>;
// a complete, view-independent description of an assembly // a complete, view-independent description of an assembly
#[derive(Clone)] #[derive(Clone)]
pub struct Assembly { pub struct Assembly {
@ -134,16 +126,7 @@ pub struct Assembly {
pub elements: Signal<Slab<Element>>, pub elements: Signal<Slab<Element>>,
pub constraints: Signal<Slab<Constraint>>, pub constraints: Signal<Slab<Constraint>>,
// solution variety tangent space. the basis vectors are stored in // solution variety tangent space
// configuration matrix format, ordered according to the elements' column
// indices. when you realize the assembly, every element that's present
// during realization gets a column index and is reflected in the tangent
// space. since the methods in this module never assign column indices
// without later realizing the assembly, we get the following invariant:
//
// (1) if an element has a column index, its tangent motions can be found
// in that column of the tangent space basis matrices
//
pub tangent: Signal<ConfigSubspace>, pub tangent: Signal<ConfigSubspace>,
// indexing // indexing
@ -155,7 +138,7 @@ impl Assembly {
Assembly { Assembly {
elements: create_signal(Slab::new()), elements: create_signal(Slab::new()),
constraints: create_signal(Slab::new()), constraints: create_signal(Slab::new()),
tangent: create_signal(ConfigSubspace::zero(0)), tangent: create_signal(ConfigSubspace::zero()),
elements_by_id: create_signal(FxHashMap::default()) elements_by_id: create_signal(FxHashMap::default())
} }
} }
@ -169,6 +152,13 @@ impl Assembly {
let id = elt.id.clone(); let id = elt.id.clone();
let key = self.elements.update(|elts| elts.insert(elt)); let key = self.elements.update(|elts| elts.insert(elt));
self.elements_by_id.update(|elts_by_id| elts_by_id.insert(id, key)); self.elements_by_id.update(|elts_by_id| elts_by_id.insert(id, key));
// realize to update the tangent space
/* KLUDGE */
// since the newly inserted element is unconstrained, we should be able
// to update the tangent space without recomputing the Hessian and its
// eigendecomposition
self.realize();
} }
pub fn try_insert_element(&self, elt: Element) -> bool { pub fn try_insert_element(&self, elt: Element) -> bool {
@ -219,7 +209,7 @@ impl Assembly {
// index the elements // index the elements
self.elements.update_silent(|elts| { self.elements.update_silent(|elts| {
for (index, (_, elt)) in elts.into_iter().enumerate() { for (index, (_, elt)) in elts.into_iter().enumerate() {
elt.column_index = Some(index); elt.column_index = index;
} }
}); });
@ -231,8 +221,8 @@ impl Assembly {
for (_, cst) in csts { for (_, cst) in csts {
if cst.active.get_untracked() && cst.lorentz_prod_valid.get_untracked() { if cst.active.get_untracked() && cst.lorentz_prod_valid.get_untracked() {
let subjects = cst.subjects; let subjects = cst.subjects;
let row = elts[subjects.0].column_index.unwrap(); let row = elts[subjects.0].column_index;
let col = elts[subjects.1].column_index.unwrap(); let col = elts[subjects.1].column_index;
gram_to_be.push_sym(row, col, cst.lorentz_prod.get_untracked()); gram_to_be.push_sym(row, col, cst.lorentz_prod.get_untracked());
} }
} }
@ -242,7 +232,7 @@ impl Assembly {
// Gram matrix // Gram matrix
let mut guess_to_be = DMatrix::<f64>::zeros(5, elts.len()); let mut guess_to_be = DMatrix::<f64>::zeros(5, elts.len());
for (_, elt) in elts { for (_, elt) in elts {
let index = elt.column_index.unwrap(); let index = elt.column_index;
gram_to_be.push_sym(index, index, 1.0); gram_to_be.push_sym(index, index, 1.0);
guess_to_be.set_column(index, &elt.representation.get_clone_untracked()); guess_to_be.set_column(index, &elt.representation.get_clone_untracked());
} }
@ -289,7 +279,7 @@ impl Assembly {
// read out the solution // read out the solution
for (_, elt) in self.elements.get_clone_untracked() { for (_, elt) in self.elements.get_clone_untracked() {
elt.representation.update( elt.representation.update(
|rep| rep.set_column(0, &config.column(elt.column_index.unwrap())) |rep| rep.set_column(0, &config.column(elt.column_index))
); );
} }
@ -300,67 +290,26 @@ impl Assembly {
// --- deformation --- // --- deformation ---
// project the given motion to the tangent space of the solution variety and pub fn deform(&self, element_motions: Vec<(ElementKey, DVectorView<f64>)>) {
// move the assembly along it. the implementation is based on invariant (1)
// from above and the following additional invariant:
//
// (2) if an element is affected by a constraint, it has a column index
//
// we have this invariant because the assembly gets realized each time you
// add a constraint
pub fn deform(&self, motion: AssemblyMotion) {
/* KLUDGE */ /* KLUDGE */
// when the tangent space is zero, deformation won't do anything, but // when the tangent space is zero, we currently need to avoid calling
// the attempt to deform should be registered in the UI. this console // its `proj` method, because it will panic rather than returning zero.
// message will do for now // in the future, we'll want a more intentionally designed system for
if self.tangent.with(|tan| tan.dim() <= 0 && tan.assembly_dim() > 0) { // handling this case
if self.tangent.with(|tan| tan.dim() <= 0) {
console::log_1(&JsValue::from("The assembly is rigid")); console::log_1(&JsValue::from("The assembly is rigid"));
return;
} }
// give a column index to each moving element that doesn't have one yet. const ELEMENT_DIM: usize = 5;
// this temporarily breaks invariant (1), but the invariant will be let assembly_dim = self.elements.with(|elts| elts.len());
// restored when we realize the assembly at the end of the deformation. let mut motion_proj = DMatrix::zeros(ELEMENT_DIM, assembly_dim);
// in the process, we find out how many matrix columns we'll need to
// hold the deformation
let realized_dim = self.tangent.with(|tan| tan.assembly_dim());
let motion_dim = self.elements.update_silent(|elts| {
let mut next_column_index = realized_dim;
for elt_motion in motion.iter() {
let moving_elt = &mut elts[elt_motion.key];
if moving_elt.column_index.is_none() {
moving_elt.column_index = Some(next_column_index);
next_column_index += 1;
}
}
next_column_index
});
// project the element motions onto the tangent space of the solution // project the element motions onto the tangent space of the solution
// variety and sum them to get a deformation of the whole assembly. the // variety, and sum them to get a deformation of the whole assembly
// matrix `motion_proj` that holds the deformation has extra columns for for (elt_key, elt_motion) in element_motions {
// any moving elements that aren't reflected in the saved tangent space let column_index = self.elements.with(|elts| elts[elt_key].column_index);
const ELEMENT_DIM: usize = 5; motion_proj += self.tangent.with(|tan| tan.proj(&elt_motion, column_index));
let mut motion_proj = DMatrix::zeros(ELEMENT_DIM, motion_dim);
for elt_motion in motion {
// we can unwrap the column index because we know that every moving
// element has one at this point
let column_index = self.elements.with_untracked(
|elts| elts[elt_motion.key].column_index.unwrap()
);
if column_index < realized_dim {
// this element had a column index when we started, so by
// invariant (1), it's reflected in the tangent space
let mut target_columns = motion_proj.columns_mut(0, realized_dim);
target_columns += self.tangent.with(
|tan| tan.proj(&elt_motion.velocity, column_index)
);
} else {
// 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;
}
} }
// step each element along the mass shell geodesic that matches its // step each element along the mass shell geodesic that matches its
@ -370,24 +319,13 @@ impl Assembly {
// element is on the 1 mass shell // element is on the 1 mass shell
for (_, elt) in self.elements.get_clone_untracked() { for (_, elt) in self.elements.get_clone_untracked() {
elt.representation.update_silent(|rep| { elt.representation.update_silent(|rep| {
match elt.column_index { let rep_next = &*rep + motion_proj.column(elt.column_index);
Some(column_index) => {
let rep_next = &*rep + motion_proj.column(column_index);
let normalizer = rep_next.dot(&(&*Q * &rep_next)); let normalizer = rep_next.dot(&(&*Q * &rep_next));
rep.set_column(0, &(rep_next / normalizer)); rep.set_column(0, &(rep_next / normalizer));
},
None => {
console::log_1(&JsValue::from(
format!("No velocity to unpack for fresh element \"{}\"", elt.id)
))
}
};
}); });
} }
// bring the configuration back onto the solution variety. this also // bring the configuration back onto the solution variety
// gets the elements' column indices and the saved tangent space back in
// sync
self.realize(); self.realize();
} }
} }

View file

@ -14,7 +14,7 @@ use web_sys::{
wasm_bindgen::{JsCast, JsValue} wasm_bindgen::{JsCast, JsValue}
}; };
use crate::{AppState, assembly::{ElementKey, ElementMotion}}; use crate::{AppState, assembly::ElementKey};
fn compile_shader( fn compile_shader(
context: &WebGl2RenderingContext, context: &WebGl2RenderingContext,
@ -341,14 +341,7 @@ pub fn Display() -> View {
]) ])
}; };
let elt_motion: DVector<f64> = time_step * TRANSLATION_SPEED * vel_field * rep; let elt_motion: DVector<f64> = time_step * TRANSLATION_SPEED * vel_field * rep;
assembly_for_raf.deform( assembly_for_raf.deform(vec![(sel, elt_motion.as_view())]);
vec![
ElementMotion {
key: sel,
velocity: elt_motion.as_view()
}
]
);
scene_changed.set(true); scene_changed.set(true);
} }
} }

View file

@ -88,17 +88,11 @@ impl PartialMatrix {
// --- configuration subspaces --- // --- configuration subspaces ---
#[derive(Clone)] #[derive(Clone)]
pub struct ConfigSubspace { pub struct ConfigSubspace(Vec<DMatrix<f64>>);
assembly_dim: usize,
basis: Vec<DMatrix<f64>>
}
impl ConfigSubspace { impl ConfigSubspace {
pub fn zero(assembly_dim: usize) -> ConfigSubspace { pub fn zero() -> ConfigSubspace {
ConfigSubspace { ConfigSubspace(Vec::new())
assembly_dim: assembly_dim,
basis: Vec::new()
}
} }
// approximate the kernel of a symmetric endomorphism of the configuration // approximate the kernel of a symmetric endomorphism of the configuration
@ -125,33 +119,26 @@ impl ConfigSubspace {
format!("Eigenvalues used to find kernel: {}", eig.eigenvalues) format!("Eigenvalues used to find kernel: {}", eig.eigenvalues)
)); ));
ConfigSubspace { ConfigSubspace(basis.collect())
assembly_dim: assembly_dim,
basis: basis.collect()
}
} }
pub fn dim(&self) -> usize { pub fn dim(&self) -> usize {
self.basis.len() let ConfigSubspace(basis) = self;
} basis.len()
pub fn assembly_dim(&self) -> usize {
self.assembly_dim
} }
// find the projection onto this subspace of the motion where the element // find the projection onto this subspace of the motion where the element
// with the given column index has velocity `v` // with the given column index has velocity `v`
/* TO DO */
// for the zero subspace, this method's behavior doesn't match its name: it
// panics rather than returning zero
pub fn proj(&self, v: &DVectorView<f64>, column_index: usize) -> DMatrix<f64> { pub fn proj(&self, v: &DVectorView<f64>, column_index: usize) -> DMatrix<f64> {
if self.dim() == 0 { let ConfigSubspace(basis) = self;
const ELEMENT_DIM: usize = 5; basis.into_iter().map(
DMatrix::zeros(ELEMENT_DIM, self.assembly_dim)
} else {
self.basis.iter().map(
|b| b.column(column_index).dot(&v) * b |b| b.column(column_index).dot(&v) * b
).sum() ).sum()
} }
} }
}
// --- descent history --- // --- descent history ---
@ -338,14 +325,14 @@ pub fn realize_gram(
state = better_state; state = better_state;
history.backoff_steps.push(backoff_steps); history.backoff_steps.push(backoff_steps);
}, },
None => return (state.config, ConfigSubspace::zero(assembly_dim), false, history) None => return (state.config, ConfigSubspace::zero(), false, history)
}; };
} }
let success = state.loss < tol; let success = state.loss < tol;
let tangent = if success { let tangent = if success {
ConfigSubspace::symmetric_kernel(hess, assembly_dim) ConfigSubspace::symmetric_kernel(hess, assembly_dim)
} else { } else {
ConfigSubspace::zero(assembly_dim) ConfigSubspace::zero()
}; };
(state.config, tangent, success, history) (state.config, tangent, success, history)
} }