Compare commits
11 commits
dde499b736
...
19e5458e1b
Author | SHA1 | Date | |
---|---|---|---|
19e5458e1b | |||
a4101ab81c | |||
6dbbe2ce2d | |||
627cea455c | |||
50c51ca08f | |||
c0e6ebf3d6 | |||
![]() |
c081f1a809 | ||
![]() |
adc60ac5c1 | ||
![]() |
27edbfb010 | ||
46ffd6c285 | |||
ecbbe2068c |
15 changed files with 804 additions and 513 deletions
64
README.md
64
README.md
|
@ -12,11 +12,11 @@ Note that currently this is just the barest beginnings of the project, more of a
|
|||
|
||||
### Implementation goals
|
||||
|
||||
* Comfortable, intuitive UI
|
||||
* Provide a comfortable, intuitive UI
|
||||
|
||||
* Able to run in browser (so implemented in WASM-compatible language)
|
||||
* Allow execution in browser (so implemented in WASM-compatible language)
|
||||
|
||||
* Produce scalable graphics of 3D diagrams, and maybe STL files (or other fabricatable file format) as well.
|
||||
* Produce scalable graphics of 3D diagrams, and maybe STL files (or other fabricatable file format) as well
|
||||
|
||||
## Prototype
|
||||
|
||||
|
@ -24,38 +24,40 @@ The latest prototype is in the folder `app-proto`. It includes both a user inter
|
|||
|
||||
### Install the prerequisites
|
||||
|
||||
1. Install [`rustup`](https://rust-lang.github.io/rustup/): the officially recommended Rust toolchain manager
|
||||
- It's available on Ubuntu as a [Snap](https://snapcraft.io/rustup)
|
||||
2. Call `rustup default stable` to "download the latest stable release of Rust and set it as your default toolchain"
|
||||
- If you forget, the `rustup` [help system](https://github.com/rust-lang/rustup/blob/d9b3601c3feb2e88cf3f8ca4f7ab4fdad71441fd/src/errors.rs#L109-L112) will remind you
|
||||
3. Call `rustup target add wasm32-unknown-unknown` to add the [most generic 32-bit WebAssembly target](https://doc.rust-lang.org/nightly/rustc/platform-support/wasm32-unknown-unknown.html)
|
||||
4. Call `cargo install wasm-pack` to install the [WebAssembly toolchain](https://rustwasm.github.io/docs/wasm-pack/)
|
||||
5. Call `cargo install trunk` to install the [Trunk](https://trunkrs.dev/) web-build tool
|
||||
1. Install [`rustup`](https://rust-lang.github.io/rustup/): the officially recommended Rust toolchain manager.
|
||||
- It's available on Ubuntu as a [Snap](https://snapcraft.io/rustup).
|
||||
2. Call `rustup default stable` to "download the latest stable release of Rust and set it as your default toolchain".
|
||||
- If you forget, the `rustup` [help system](https://github.com/rust-lang/rustup/blob/d9b3601c3feb2e88cf3f8ca4f7ab4fdad71441fd/src/errors.rs#L109-L112) will remind you.
|
||||
3. Call `rustup target add wasm32-unknown-unknown` to add the [most generic 32-bit WebAssembly target](https://doc.rust-lang.org/nightly/rustc/platform-support/wasm32-unknown-unknown.html).
|
||||
4. Call `cargo install wasm-pack` to install the [WebAssembly toolchain](https://rustwasm.github.io/docs/wasm-pack/).
|
||||
5. Call `cargo install trunk` to install the [Trunk](https://trunkrs.dev/) web-build tool.
|
||||
- In the future, `trunk` can be updated with the same command. (You may need the `--locked` flag if your ambient version of `rustc` does not match that required by `trunk`.)
|
||||
6. Add the `.cargo/bin` folder in your home directory to your executable search path
|
||||
- This lets you call Trunk, and other tools installed by Cargo, without specifying their paths
|
||||
- On POSIX systems, the search path is stored in the `PATH` environment variable
|
||||
- This lets you call Trunk, and other tools installed by Cargo, without specifying their paths.
|
||||
- On POSIX systems, the search path is stored in the `PATH` environment variable.
|
||||
- Alternatively, if you don't want to adjust your `PATH`, you can install `trunk` in another directory `DIR` via `cargo install --root DIR trunk`.
|
||||
|
||||
### Play with the prototype
|
||||
|
||||
1. From the `app-proto` folder, call `trunk serve --release` to build and serve the prototype
|
||||
- The crates the prototype depends on will be downloaded and served automatically
|
||||
- For a faster build, at the expense of a much slower prototype, you can call `trunk serve` without the `--release` flag
|
||||
1. From the `app-proto` folder, call `trunk serve --release` to build and serve the prototype.
|
||||
- The crates the prototype depends on will be downloaded and served automatically.
|
||||
- For a faster build, at the expense of a much slower prototype, you can call `trunk serve` without the `--release` flag.
|
||||
- If you want to stay in the top-level folder, you can call `trunk serve --config app-proto [--release]` from there instead.
|
||||
3. In a web browser, visit one of the URLs listed under the message `INFO 📡 server listening at:`
|
||||
- Touching any file in the `app-proto` folder will make Trunk rebuild and live-reload the prototype
|
||||
4. Press *ctrl+C* in the shell where Trunk is running to stop serving the prototype
|
||||
3. In a web browser, visit one of the URLs listed under the message `INFO 📡 server listening at:`.
|
||||
- Touching any file in the `app-proto` folder will make Trunk rebuild and live-reload the prototype.
|
||||
4. Press *ctrl+C* in the shell where Trunk is running to stop serving the prototype.
|
||||
|
||||
### Run the engine on some example problems
|
||||
|
||||
1. Use `sh` to run the script `tools/run-examples.sh`
|
||||
- The script is location-independent, so you can do this from anywhere in the dyna3 repository
|
||||
1. Use `sh` to run the script `tools/run-examples.sh`.
|
||||
- The script is location-independent, so you can do this from anywhere in the dyna3 repository.
|
||||
- The call from the top level of the repository is:
|
||||
|
||||
```bash
|
||||
sh tools/run-examples.sh
|
||||
```
|
||||
- For each example problem, the engine will print the value of the loss function at each optimization step
|
||||
- The first example that prints is the same as the Irisawa hexlet example from the Julia version of the engine prototype. If you go into `engine-proto/gram-test`, launch Julia, and then
|
||||
- For each example problem, the engine will print the value of the loss function at each optimization step.
|
||||
- The first example that prints is the same as the Irisawa hexlet example from the Julia version of the engine prototype. If you go into `engine-proto/gram-test`, launch Julia, and then execute
|
||||
|
||||
```julia
|
||||
include("irisawa-hexlet.jl")
|
||||
|
@ -64,24 +66,24 @@ The latest prototype is in the folder `app-proto`. It includes both a user inter
|
|||
end
|
||||
```
|
||||
|
||||
you should see that it prints basically the same loss history until the last few steps, when the lower default precision of the Rust engine really starts to show
|
||||
you should see that it prints basically the same loss history until the last few steps, when the lower default precision of the Rust engine really starts to show.
|
||||
|
||||
### Run the automated tests
|
||||
|
||||
1. Go into the `app-proto` folder
|
||||
2. Call `cargo test`
|
||||
1. Go into the `app-proto` folder.
|
||||
2. Call `cargo test`.
|
||||
|
||||
### Deploy the prototype
|
||||
|
||||
1. From the `app-proto` folder, call `trunk build --release`
|
||||
- Building in [release mode](https://doc.rust-lang.org/cargo/reference/profiles.html#release) produces an executable which is smaller and often much faster, but harder to debug and more time-consuming to build
|
||||
- If you want to stay in the top-level folder, you can call `trunk build --config app-proto --release` from there instead
|
||||
1. From the `app-proto` folder, call `trunk build --release`.
|
||||
- Building in [release mode](https://doc.rust-lang.org/cargo/reference/profiles.html#release) produces an executable which is smaller and often much faster, but harder to debug and more time-consuming to build.
|
||||
- If you want to stay in the top-level folder, you can call `trunk build --config app-proto --release` from there instead.
|
||||
2. Use `sh` to run the packaging script `tools/package-for-deployment.sh`.
|
||||
- The script is location-independent, so you can do this from anywhere in the dyna3 repository
|
||||
- The script is location-independent, so you can do this from anywhere in the dyna3 repository.
|
||||
- The call from the top level of the repository is:
|
||||
```bash
|
||||
sh tools/package-for-deployment.sh
|
||||
```
|
||||
- This will overwrite or replace the files in `deploy/dyna3`
|
||||
- This will overwrite or replace the files in `deploy/dyna3`.
|
||||
3. Put the contents of `deploy/dyna3` in the folder on your server that the prototype will be served from.
|
||||
- To simplify uploading, you might want to combine these files into an archive called `deploy/dyna3.zip`. Git has been set to ignore this path
|
||||
- To simplify uploading, you might want to combine these files into an archive called `deploy/dyna3.zip`. Git has been set to ignore this path.
|
||||
|
|
21
app-proto/Cargo.lock
generated
21
app-proto/Cargo.lock
generated
|
@ -255,6 +255,7 @@ dependencies = [
|
|||
"charming",
|
||||
"console_error_panic_hook",
|
||||
"dyna3",
|
||||
"enum-iterator",
|
||||
"itertools",
|
||||
"js-sys",
|
||||
"lazy_static",
|
||||
|
@ -271,6 +272,26 @@ version = "1.13.0"
|
|||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "60b1af1c220855b6ceac025d3f6ecdd2b7c4894bfe9cd9bda4fbb4bc7c0d4cf0"
|
||||
|
||||
[[package]]
|
||||
name = "enum-iterator"
|
||||
version = "2.3.0"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "a4549325971814bda7a44061bf3fe7e487d447cba01e4220a4b454d630d7a016"
|
||||
dependencies = [
|
||||
"enum-iterator-derive",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "enum-iterator-derive"
|
||||
version = "1.5.0"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "685adfa4d6f3d765a26bc5dbc936577de9abf756c1feeb3089b01dd395034842"
|
||||
dependencies = [
|
||||
"proc-macro2",
|
||||
"quote",
|
||||
"syn",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "equivalent"
|
||||
version = "1.0.1"
|
||||
|
|
|
@ -10,6 +10,7 @@ default = ["console_error_panic_hook"]
|
|||
dev = []
|
||||
|
||||
[dependencies]
|
||||
enum-iterator = "2.3.0"
|
||||
itertools = "0.13.0"
|
||||
js-sys = "0.3.70"
|
||||
lazy_static = "1.5.0"
|
||||
|
|
|
@ -1,10 +1,11 @@
|
|||
use enum_iterator::{all, Sequence};
|
||||
use nalgebra::{DMatrix, DVector, DVectorView};
|
||||
use std::{
|
||||
cell::Cell,
|
||||
cmp::Ordering,
|
||||
collections::{BTreeMap, BTreeSet},
|
||||
fmt,
|
||||
fmt::{Debug, Formatter},
|
||||
fmt::{Debug, Display, Formatter},
|
||||
hash::{Hash, Hasher},
|
||||
rc::Rc,
|
||||
sync::{atomic, atomic::AtomicU64},
|
||||
|
@ -26,6 +27,7 @@ use crate::{
|
|||
ConfigSubspace,
|
||||
ConstraintProblem,
|
||||
DescentHistory,
|
||||
MatrixEntry,
|
||||
Realization,
|
||||
},
|
||||
specified::SpecifiedValue,
|
||||
|
@ -84,6 +86,14 @@ impl Ord for dyn Serial {
|
|||
}
|
||||
}
|
||||
|
||||
// Small helper function to generate consistent errors when there
|
||||
// are indexing issues in a ProblemPoser
|
||||
fn indexing_error(item: &str, name: &str, actor: &str) -> String {
|
||||
format!(
|
||||
"{item} \"{name}\" must be indexed before {actor} writes problem data"
|
||||
)
|
||||
}
|
||||
|
||||
pub trait ProblemPoser {
|
||||
fn pose(&self, problem: &mut ConstraintProblem);
|
||||
}
|
||||
|
@ -125,8 +135,8 @@ pub trait Element: Serial + ProblemPoser + DisplayItem {
|
|||
}
|
||||
|
||||
impl Debug for dyn Element {
|
||||
fn fmt(&self, f: &mut Formatter<'_>) -> Result<(), fmt::Error> {
|
||||
self.id().fmt(f)
|
||||
fn fmt(&self, f: &mut Formatter<'_>) -> fmt::Result {
|
||||
Debug::fmt(&self.id(), f)
|
||||
}
|
||||
}
|
||||
|
||||
|
@ -249,10 +259,10 @@ impl Serial for Sphere {
|
|||
impl ProblemPoser for Sphere {
|
||||
fn pose(&self, problem: &mut ConstraintProblem) {
|
||||
let index = self.column_index().expect(
|
||||
format!("Sphere \"{}\" should be indexed before writing problem data", self.id).as_str()
|
||||
);
|
||||
indexing_error("Sphere", &self.id, "it").as_str());
|
||||
problem.gram.push_sym(index, index, 1.0);
|
||||
problem.guess.set_column(index, &self.representation.get_clone_untracked());
|
||||
problem.guess.set_column(
|
||||
index, &self.representation.get_clone_untracked());
|
||||
}
|
||||
}
|
||||
|
||||
|
@ -269,6 +279,7 @@ pub struct Point {
|
|||
|
||||
impl Point {
|
||||
const WEIGHT_COMPONENT: usize = 3;
|
||||
const NORM_COMPONENT: usize = 4;
|
||||
|
||||
pub fn new(
|
||||
id: String,
|
||||
|
@ -303,6 +314,15 @@ impl Element for Point {
|
|||
)
|
||||
}
|
||||
|
||||
fn default_regulators(self: Rc<Self>) -> Vec<Rc<dyn Regulator>> {
|
||||
all::<Axis>()
|
||||
.map(|axis| {
|
||||
Rc::new(PointCoordinateRegulator::new(self.clone(), axis))
|
||||
as Rc::<dyn Regulator>
|
||||
})
|
||||
.collect()
|
||||
}
|
||||
|
||||
fn id(&self) -> &String {
|
||||
&self.id
|
||||
}
|
||||
|
@ -345,11 +365,11 @@ impl Serial for Point {
|
|||
impl ProblemPoser for Point {
|
||||
fn pose(&self, problem: &mut ConstraintProblem) {
|
||||
let index = self.column_index().expect(
|
||||
format!("Point \"{}\" should be indexed before writing problem data", self.id).as_str()
|
||||
);
|
||||
indexing_error("Point", &self.id, "it").as_str());
|
||||
problem.gram.push_sym(index, index, 0.0);
|
||||
problem.frozen.push(Self::WEIGHT_COMPONENT, index, 0.5);
|
||||
problem.guess.set_column(index, &self.representation.get_clone_untracked());
|
||||
problem.guess.set_column(
|
||||
index, &self.representation.get_clone_untracked());
|
||||
}
|
||||
}
|
||||
|
||||
|
@ -394,7 +414,8 @@ pub struct InversiveDistanceRegulator {
|
|||
|
||||
impl InversiveDistanceRegulator {
|
||||
pub fn new(subjects: [Rc<dyn Element>; 2]) -> Self {
|
||||
let representations = subjects.each_ref().map(|subj| subj.representation());
|
||||
let representations = subjects.each_ref().map(
|
||||
|subj| subj.representation());
|
||||
let measurement = create_memo(move || {
|
||||
representations[0].with(|rep_0|
|
||||
representations[1].with(|rep_1|
|
||||
|
@ -436,8 +457,8 @@ impl ProblemPoser for InversiveDistanceRegulator {
|
|||
if let Some(val) = set_pt.value {
|
||||
let [row, col] = self.subjects.each_ref().map(
|
||||
|subj| subj.column_index().expect(
|
||||
"Subjects should be indexed before inversive distance regulator writes problem data"
|
||||
)
|
||||
indexing_error("Subject", subj.id(),
|
||||
"inversive distance regulator").as_str())
|
||||
);
|
||||
problem.gram.push_sym(row, col, val);
|
||||
}
|
||||
|
@ -446,14 +467,14 @@ impl ProblemPoser for InversiveDistanceRegulator {
|
|||
}
|
||||
|
||||
pub struct HalfCurvatureRegulator {
|
||||
pub subject: Rc<dyn Element>,
|
||||
pub subject: Rc<Sphere>,
|
||||
pub measurement: ReadSignal<f64>,
|
||||
pub set_point: Signal<SpecifiedValue>,
|
||||
serial: u64,
|
||||
}
|
||||
|
||||
impl HalfCurvatureRegulator {
|
||||
pub fn new(subject: Rc<dyn Element>) -> Self {
|
||||
pub fn new(subject: Rc<Sphere>) -> Self {
|
||||
let measurement = subject.representation().map(
|
||||
|rep| rep[Sphere::CURVATURE_COMPONENT]
|
||||
);
|
||||
|
@ -490,14 +511,87 @@ impl ProblemPoser for HalfCurvatureRegulator {
|
|||
self.set_point.with_untracked(|set_pt| {
|
||||
if let Some(val) = set_pt.value {
|
||||
let col = self.subject.column_index().expect(
|
||||
"Subject should be indexed before half-curvature regulator writes problem data"
|
||||
);
|
||||
indexing_error("Subject", &self.subject.id,
|
||||
"half-curvature regulator").as_str());
|
||||
problem.frozen.push(Sphere::CURVATURE_COMPONENT, col, val);
|
||||
}
|
||||
});
|
||||
}
|
||||
}
|
||||
|
||||
#[derive(Clone, Copy, Sequence)]
|
||||
pub enum Axis { X = 0, Y = 1, Z = 2 }
|
||||
|
||||
impl Axis {
|
||||
fn name(&self) -> &'static str {
|
||||
match self { Axis::X => "X", Axis::Y => "Y", Axis::Z => "Z" }
|
||||
}
|
||||
}
|
||||
|
||||
impl Display for Axis {
|
||||
fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
|
||||
f.write_str(self.name())
|
||||
}
|
||||
}
|
||||
|
||||
pub struct PointCoordinateRegulator {
|
||||
pub subject: Rc<Point>,
|
||||
pub axis: Axis,
|
||||
pub measurement: ReadSignal<f64>,
|
||||
pub set_point: Signal<SpecifiedValue>,
|
||||
serial: u64
|
||||
}
|
||||
|
||||
impl PointCoordinateRegulator {
|
||||
pub fn new(subject: Rc<Point>, axis: Axis) -> Self {
|
||||
let measurement = subject.representation().map(
|
||||
move |rep| rep[axis as usize]
|
||||
);
|
||||
let set_point = create_signal(SpecifiedValue::from_empty_spec());
|
||||
Self {
|
||||
subject, axis, measurement, set_point, serial: Self::next_serial()
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
impl Serial for PointCoordinateRegulator {
|
||||
fn serial(&self) -> u64 { self.serial }
|
||||
}
|
||||
|
||||
impl Regulator for PointCoordinateRegulator {
|
||||
fn subjects(&self) -> Vec<Rc<dyn Element>> { vec![self.subject.clone()] }
|
||||
fn measurement(&self) -> ReadSignal<f64> { self.measurement }
|
||||
fn set_point(&self) -> Signal<SpecifiedValue> { self.set_point }
|
||||
}
|
||||
|
||||
impl ProblemPoser for PointCoordinateRegulator {
|
||||
fn pose(&self, problem: &mut ConstraintProblem) {
|
||||
self.set_point.with_untracked(|set_pt| {
|
||||
if let Some(val) = set_pt.value {
|
||||
let col = self.subject.column_index().expect(
|
||||
indexing_error("Subject", &self.subject.id,
|
||||
"point-coordinate regulator").as_str());
|
||||
problem.frozen.push(self.axis as usize, col, val);
|
||||
// If all three of the subject's spatial coordinates have been
|
||||
// frozen, then freeze its norm component:
|
||||
let mut coords = [0.0; Axis::CARDINALITY];
|
||||
let mut nset: usize = 0;
|
||||
for &MatrixEntry {index, value} in &(problem.frozen) {
|
||||
if index.1 == col && index.0 < Axis::CARDINALITY {
|
||||
nset += 1;
|
||||
coords[index.0] = value
|
||||
}
|
||||
}
|
||||
if nset == Axis::CARDINALITY {
|
||||
let [x, y, z] = coords;
|
||||
problem.frozen.push(Point::NORM_COMPONENT,
|
||||
col, point(x,y,z)[Point::NORM_COMPONENT]);
|
||||
}
|
||||
}
|
||||
});
|
||||
}
|
||||
}
|
||||
|
||||
// the velocity is expressed in uniform coordinates
|
||||
pub struct ElementMotion<'a> {
|
||||
pub element: Rc<dyn Element>,
|
||||
|
@ -590,7 +684,8 @@ impl Assembly {
|
|||
let id = elt.id().clone();
|
||||
let elt_rc = Rc::new(elt);
|
||||
self.elements.update(|elts| elts.insert(elt_rc.clone()));
|
||||
self.elements_by_id.update(|elts_by_id| elts_by_id.insert(id, elt_rc.clone()));
|
||||
self.elements_by_id.update(
|
||||
|elts_by_id| elts_by_id.insert(id, elt_rc.clone()));
|
||||
|
||||
// create and insert the element's default regulators
|
||||
for reg in elt_rc.default_regulators() {
|
||||
|
@ -666,7 +761,8 @@ impl Assembly {
|
|||
pub fn load_config(&self, config: &DMatrix<f64>) {
|
||||
for elt in self.elements.get_clone_untracked() {
|
||||
elt.representation().update(
|
||||
|rep| rep.set_column(0, &config.column(elt.column_index().unwrap()))
|
||||
|rep| rep.set_column(
|
||||
0, &config.column(elt.column_index().unwrap()))
|
||||
);
|
||||
}
|
||||
}
|
||||
|
@ -698,6 +794,7 @@ impl Assembly {
|
|||
/* DEBUG */
|
||||
// log the Gram matrix
|
||||
console_log!("Gram matrix:\n{}", problem.gram);
|
||||
console_log!("Frozen entries:\n{}", problem.frozen);
|
||||
|
||||
/* DEBUG */
|
||||
// log the initial configuration matrix
|
||||
|
@ -810,7 +907,8 @@ impl Assembly {
|
|||
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);
|
||||
let mut target_columns =
|
||||
motion_proj.columns_mut(0, realized_dim);
|
||||
target_columns += self.tangent.with(
|
||||
|tan| tan.proj(&elt_motion.velocity, column_index)
|
||||
);
|
||||
|
@ -818,7 +916,8 @@ 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);
|
||||
let unif_to_std = elt_motion.element.representation().with_untracked(
|
||||
let unif_to_std =
|
||||
elt_motion.element.representation().with_untracked(
|
||||
|rep| local_unif_to_std(rep.as_view())
|
||||
);
|
||||
target_column += unif_to_std * elt_motion.velocity;
|
||||
|
@ -837,7 +936,10 @@ impl Assembly {
|
|||
elt.project_to_normalized(rep);
|
||||
},
|
||||
None => {
|
||||
console_log!("No velocity to unpack for fresh element \"{}\"", elt.id())
|
||||
console_log!(
|
||||
"No velocity to unpack for fresh element \"{}\"",
|
||||
elt.id()
|
||||
)
|
||||
},
|
||||
};
|
||||
});
|
||||
|
@ -857,7 +959,8 @@ mod tests {
|
|||
use crate::engine;
|
||||
|
||||
#[test]
|
||||
#[should_panic(expected = "Sphere \"sphere\" should be indexed before writing problem data")]
|
||||
#[should_panic(expected =
|
||||
"Sphere \"sphere\" must be indexed before it writes problem data")]
|
||||
fn unindexed_element_test() {
|
||||
let _ = create_root(|| {
|
||||
let elt = Sphere::default("sphere".to_string(), 0);
|
||||
|
@ -866,17 +969,20 @@ mod tests {
|
|||
}
|
||||
|
||||
#[test]
|
||||
#[should_panic(expected = "Subjects should be indexed before inversive distance regulator writes problem data")]
|
||||
#[should_panic(expected = "Subject \"sphere1\" must be indexed before \
|
||||
inversive distance regulator writes problem data")]
|
||||
fn unindexed_subject_test_inversive_distance() {
|
||||
let _ = create_root(|| {
|
||||
let subjects = [0, 1].map(
|
||||
|k| Rc::new(Sphere::default(format!("sphere{k}"), k)) as Rc<dyn Element>
|
||||
|k| Rc::new(
|
||||
Sphere::default(format!("sphere{k}"), k)) as Rc<dyn Element>
|
||||
);
|
||||
subjects[0].set_column_index(0);
|
||||
InversiveDistanceRegulator {
|
||||
subjects: subjects,
|
||||
measurement: create_memo(|| 0.0),
|
||||
set_point: create_signal(SpecifiedValue::try_from("0.0".to_string()).unwrap()),
|
||||
set_point: create_signal(
|
||||
SpecifiedValue::try_from("0.0".to_string()).unwrap()),
|
||||
serial: InversiveDistanceRegulator::next_serial()
|
||||
}.pose(&mut ConstraintProblem::new(2));
|
||||
});
|
||||
|
@ -905,8 +1011,10 @@ mod tests {
|
|||
// nudge the sphere repeatedly along the `z` axis
|
||||
const STEP_SIZE: f64 = 0.0025;
|
||||
const STEP_CNT: usize = 400;
|
||||
let sphere = assembly.elements_by_id.with(|elts_by_id| elts_by_id[sphere_id].clone());
|
||||
let velocity = DVector::from_column_slice(&[0.0, 0.0, STEP_SIZE, 0.0]);
|
||||
let sphere = assembly.elements_by_id.with(
|
||||
|elts_by_id| elts_by_id[sphere_id].clone());
|
||||
let velocity =
|
||||
DVector::from_column_slice(&[0.0, 0.0, STEP_SIZE, 0.0]);
|
||||
for _ in 0..STEP_CNT {
|
||||
assembly.deform(
|
||||
vec![
|
||||
|
@ -924,7 +1032,8 @@ mod tests {
|
|||
let final_half_curv = sphere.representation().with_untracked(
|
||||
|rep| rep[Sphere::CURVATURE_COMPONENT]
|
||||
);
|
||||
assert!((final_half_curv / INITIAL_HALF_CURV - 1.0).abs() < DRIFT_TOL);
|
||||
assert!((final_half_curv / INITIAL_HALF_CURV - 1.0).abs()
|
||||
< DRIFT_TOL);
|
||||
});
|
||||
}
|
||||
}
|
|
@ -39,7 +39,9 @@ pub fn AddRemove() -> View {
|
|||
}
|
||||
) { "Add point" }
|
||||
button(
|
||||
class = "emoji", /* KLUDGE */ // for convenience, we're using an emoji as a temporary icon for this button
|
||||
/* KLUDGE */ // for convenience, we're using an emoji as an
|
||||
// icon for this button
|
||||
class = "emoji",
|
||||
disabled = {
|
||||
let state = use_context::<AppState>();
|
||||
state.selection.with(|sel| sel.len() != 2)
|
||||
|
|
|
@ -50,7 +50,8 @@ impl SceneSpheres {
|
|||
}
|
||||
|
||||
fn len_i32(&self) -> i32 {
|
||||
self.representations.len().try_into().expect("Number of spheres must fit in a 32-bit integer")
|
||||
self.representations.len().try_into().expect(
|
||||
"Number of spheres must fit in a 32-bit integer")
|
||||
}
|
||||
|
||||
fn push(
|
||||
|
@ -127,8 +128,12 @@ impl DisplayItem for Sphere {
|
|||
const HIGHLIGHT: f32 = 0.2;
|
||||
|
||||
let representation = self.representation.get_clone_untracked();
|
||||
let color = if selected { self.color.map(|channel| 0.2 + 0.8*channel) } else { self.color };
|
||||
let opacity = if self.ghost.get() { GHOST_OPACITY } else { DEFAULT_OPACITY };
|
||||
let color =
|
||||
if selected { self.color.map(|channel| 0.2 + 0.8*channel) }
|
||||
else { self.color };
|
||||
let opacity =
|
||||
if self.ghost.get() { GHOST_OPACITY }
|
||||
else { DEFAULT_OPACITY };
|
||||
let highlight = if selected { 1.0 } else { HIGHLIGHT };
|
||||
scene.spheres.push(representation, color, opacity, highlight);
|
||||
}
|
||||
|
@ -145,7 +150,8 @@ impl DisplayItem for Sphere {
|
|||
// `a*u^2 + b*u + c` by the linear function `b*u + c`
|
||||
const DEG_THRESHOLD: f64 = 1e-9;
|
||||
|
||||
let rep = self.representation.with_untracked(|rep| assembly_to_world * rep);
|
||||
let rep = self.representation.with_untracked(
|
||||
|rep| assembly_to_world * rep);
|
||||
let a = -rep[3] * dir.norm_squared();
|
||||
let b = rep.rows_range(..3).dot(&dir);
|
||||
let c = -rep[4];
|
||||
|
@ -186,7 +192,9 @@ impl DisplayItem for Point {
|
|||
const HIGHLIGHT: f32 = 0.5;
|
||||
|
||||
let representation = self.representation.get_clone_untracked();
|
||||
let color = if selected { self.color.map(|channel| 0.2 + 0.8*channel) } else { self.color };
|
||||
let color =
|
||||
if selected { self.color.map(|channel| 0.2 + 0.8*channel) }
|
||||
else { self.color };
|
||||
let opacity = if self.ghost.get() { GHOST_OPACITY } else { 1.0 };
|
||||
let highlight = if selected { 1.0 } else { HIGHLIGHT };
|
||||
scene.points.push(representation, color, opacity, highlight, selected);
|
||||
|
@ -199,7 +207,8 @@ impl DisplayItem for Point {
|
|||
assembly_to_world: &DMatrix<f64>,
|
||||
pixel_size: f64,
|
||||
) -> Option<f64> {
|
||||
let rep = self.representation.with_untracked(|rep| assembly_to_world * rep);
|
||||
let rep = self.representation.with_untracked(
|
||||
|rep| assembly_to_world * rep);
|
||||
if rep[2] < 0.0 {
|
||||
// this constant should be kept synchronized with `point.frag`
|
||||
const POINT_RADIUS_PX: f64 = 4.0;
|
||||
|
@ -357,11 +366,12 @@ fn event_dir(event: &MouseEvent) -> (Vector3<f64>, f64) {
|
|||
// this constant should be kept synchronized with `spheres.frag` and
|
||||
// `point.vert`
|
||||
const FOCAL_SLOPE: f64 = 0.3;
|
||||
|
||||
let horizontal = f64::from(event.client_x()) - rect.left();
|
||||
let vertical = rect.bottom() - f64::from(event.client_y());
|
||||
(
|
||||
Vector3::new(
|
||||
FOCAL_SLOPE * (2.0*(f64::from(event.client_x()) - rect.left()) - width) / shortdim,
|
||||
FOCAL_SLOPE * (2.0*(rect.bottom() - f64::from(event.client_y())) - height) / shortdim,
|
||||
FOCAL_SLOPE * (2.0*horizontal - width) / shortdim,
|
||||
FOCAL_SLOPE * (2.0*vertical - height) / shortdim,
|
||||
-1.0,
|
||||
),
|
||||
FOCAL_SLOPE * 2.0 / shortdim,
|
||||
|
@ -445,7 +455,8 @@ pub fn Display() -> View {
|
|||
let performance = window().unwrap().performance().unwrap();
|
||||
|
||||
// get the display canvas
|
||||
let canvas = display.get().unchecked_into::<web_sys::HtmlCanvasElement>();
|
||||
let canvas =
|
||||
display.get().unchecked_into::<web_sys::HtmlCanvasElement>();
|
||||
let ctx = canvas
|
||||
.get_context("webgl2")
|
||||
.unwrap()
|
||||
|
@ -458,7 +469,8 @@ pub fn Display() -> View {
|
|||
|
||||
// set blend mode
|
||||
ctx.enable(WebGl2RenderingContext::BLEND);
|
||||
ctx.blend_func(WebGl2RenderingContext::SRC_ALPHA, WebGl2RenderingContext::ONE_MINUS_SRC_ALPHA);
|
||||
ctx.blend_func(WebGl2RenderingContext::SRC_ALPHA,
|
||||
WebGl2RenderingContext::ONE_MINUS_SRC_ALPHA);
|
||||
|
||||
// set up the sphere rendering program
|
||||
let sphere_program = set_up_program(
|
||||
|
@ -487,16 +499,20 @@ pub fn Display() -> View {
|
|||
// machine, the the length of a float or genType array seems to be
|
||||
// capped at 1024 elements
|
||||
console::log_2(
|
||||
&ctx.get_parameter(WebGl2RenderingContext::MAX_FRAGMENT_UNIFORM_VECTORS).unwrap(),
|
||||
&ctx.get_parameter(
|
||||
WebGl2RenderingContext::MAX_FRAGMENT_UNIFORM_VECTORS).unwrap(),
|
||||
&JsValue::from("uniform vectors available"),
|
||||
);
|
||||
|
||||
// find the sphere program's vertex attribute
|
||||
let viewport_position_attr = ctx.get_attrib_location(&sphere_program, "position") as u32;
|
||||
let viewport_position_attr =
|
||||
ctx.get_attrib_location(&sphere_program, "position") as u32;
|
||||
|
||||
// find the sphere program's uniforms
|
||||
const SPHERE_MAX: usize = 200;
|
||||
let sphere_cnt_loc = ctx.get_uniform_location(&sphere_program, "sphere_cnt");
|
||||
let sphere_cnt_loc = ctx.get_uniform_location(
|
||||
&sphere_program, "sphere_cnt"
|
||||
);
|
||||
let sphere_sp_locs = get_uniform_array_locations::<SPHERE_MAX>(
|
||||
&ctx, &sphere_program, "sphere_list", Some("sp")
|
||||
);
|
||||
|
@ -509,10 +525,18 @@ pub fn Display() -> View {
|
|||
let sphere_highlight_locs = get_uniform_array_locations::<SPHERE_MAX>(
|
||||
&ctx, &sphere_program, "highlight_list", None
|
||||
);
|
||||
let resolution_loc = ctx.get_uniform_location(&sphere_program, "resolution");
|
||||
let shortdim_loc = ctx.get_uniform_location(&sphere_program, "shortdim");
|
||||
let layer_threshold_loc = ctx.get_uniform_location(&sphere_program, "layer_threshold");
|
||||
let debug_mode_loc = ctx.get_uniform_location(&sphere_program, "debug_mode");
|
||||
let resolution_loc = ctx.get_uniform_location(
|
||||
&sphere_program, "resolution"
|
||||
);
|
||||
let shortdim_loc = ctx.get_uniform_location(
|
||||
&sphere_program, "shortdim"
|
||||
);
|
||||
let layer_threshold_loc = ctx.get_uniform_location(
|
||||
&sphere_program, "layer_threshold"
|
||||
);
|
||||
let debug_mode_loc = ctx.get_uniform_location(
|
||||
&sphere_program, "debug_mode"
|
||||
);
|
||||
|
||||
// load the viewport vertex positions into a new vertex buffer object
|
||||
const VERTEX_CNT: usize = 6;
|
||||
|
@ -526,13 +550,18 @@ pub fn Display() -> View {
|
|||
1.0, 1.0, 0.0,
|
||||
1.0, -1.0, 0.0,
|
||||
];
|
||||
let viewport_position_buffer = load_new_buffer(&ctx, &viewport_positions);
|
||||
let viewport_position_buffer =
|
||||
load_new_buffer(&ctx, &viewport_positions);
|
||||
|
||||
// find the point program's vertex attributes
|
||||
let point_position_attr = ctx.get_attrib_location(&point_program, "position") as u32;
|
||||
let point_color_attr = ctx.get_attrib_location(&point_program, "color") as u32;
|
||||
let point_highlight_attr = ctx.get_attrib_location(&point_program, "highlight") as u32;
|
||||
let point_selection_attr = ctx.get_attrib_location(&point_program, "selected") as u32;
|
||||
let point_position_attr =
|
||||
ctx.get_attrib_location(&point_program, "position") as u32;
|
||||
let point_color_attr =
|
||||
ctx.get_attrib_location(&point_program, "color") as u32;
|
||||
let point_highlight_attr =
|
||||
ctx.get_attrib_location(&point_program, "highlight") as u32;
|
||||
let point_selection_attr =
|
||||
ctx.get_attrib_location(&point_program, "selected") as u32;
|
||||
|
||||
// set up a repainting routine
|
||||
let (_, start_animation_loop, _) = create_raf(move || {
|
||||
|
@ -596,7 +625,8 @@ pub fn Display() -> View {
|
|||
let realization_successful = state.assembly.realization_status.with(
|
||||
|status| status.is_ok()
|
||||
);
|
||||
let step_val = state.assembly.step.with_untracked(|step| step.value);
|
||||
let step_val =
|
||||
state.assembly.step.with_untracked(|step| step.value);
|
||||
let on_init_step = step_val.is_some_and(|n| n == 0.0);
|
||||
let on_last_step = step_val.is_some_and(
|
||||
|n| state.assembly.descent_history.with_untracked(
|
||||
|
@ -606,7 +636,8 @@ pub fn Display() -> View {
|
|||
let on_manipulable_step =
|
||||
!realization_successful && on_init_step
|
||||
|| realization_successful && on_last_step;
|
||||
if on_manipulable_step && state.selection.with(|sel| sel.len() == 1) {
|
||||
if on_manipulable_step
|
||||
&& state.selection.with(|sel| sel.len() == 1) {
|
||||
let sel = state.selection.with(
|
||||
|sel| sel.into_iter().next().unwrap().clone()
|
||||
);
|
||||
|
@ -651,7 +682,8 @@ pub fn Display() -> View {
|
|||
// measure mean frame interval
|
||||
frames_since_last_sample += 1;
|
||||
if frames_since_last_sample >= SAMPLE_PERIOD {
|
||||
mean_frame_interval.set((time - last_sample_time) / (SAMPLE_PERIOD as f64));
|
||||
mean_frame_interval.set(
|
||||
(time - last_sample_time) / (SAMPLE_PERIOD as f64));
|
||||
last_sample_time = time;
|
||||
frames_since_last_sample = 0;
|
||||
}
|
||||
|
@ -676,7 +708,8 @@ pub fn Display() -> View {
|
|||
// set up the scene
|
||||
state.assembly.elements.with_untracked(
|
||||
|elts| for elt in elts {
|
||||
let selected = state.selection.with(|sel| sel.contains(elt));
|
||||
let selected =
|
||||
state.selection.with(|sel| sel.contains(elt));
|
||||
elt.show(&mut scene, selected);
|
||||
}
|
||||
);
|
||||
|
@ -691,7 +724,8 @@ pub fn Display() -> View {
|
|||
ctx.enable_vertex_attrib_array(viewport_position_attr);
|
||||
|
||||
// write the spheres in world coordinates
|
||||
let sphere_reps_world: Vec<_> = scene.spheres.representations.into_iter().map(
|
||||
let sphere_reps_world: Vec<_> =
|
||||
scene.spheres.representations.into_iter().map(
|
||||
|rep| (&asm_to_world * rep).cast::<f32>()
|
||||
).collect();
|
||||
|
||||
|
@ -729,10 +763,12 @@ pub fn Display() -> View {
|
|||
|
||||
// bind the viewport vertex position buffer to the position
|
||||
// attribute in the vertex shader
|
||||
bind_to_attribute(&ctx, viewport_position_attr, SPACE_DIM as i32, &viewport_position_buffer);
|
||||
bind_to_attribute(&ctx, viewport_position_attr,
|
||||
SPACE_DIM as i32, &viewport_position_buffer);
|
||||
|
||||
// draw the scene
|
||||
ctx.draw_arrays(WebGl2RenderingContext::TRIANGLES, 0, VERTEX_CNT as i32);
|
||||
ctx.draw_arrays(WebGl2RenderingContext::TRIANGLES, 0,
|
||||
VERTEX_CNT as i32);
|
||||
|
||||
// disable the sphere program's vertex attribute
|
||||
ctx.disable_vertex_attrib_array(viewport_position_attr);
|
||||
|
@ -760,13 +796,19 @@ pub fn Display() -> View {
|
|||
// load the point positions and colors into new buffers and
|
||||
// bind them to the corresponding attributes in the vertex
|
||||
// shader
|
||||
bind_new_buffer_to_attribute(&ctx, point_position_attr, SPACE_DIM as i32, point_positions.as_slice());
|
||||
bind_new_buffer_to_attribute(&ctx, point_color_attr, (COLOR_SIZE + 1) as i32, scene.points.colors_with_opacity.concat().as_slice());
|
||||
bind_new_buffer_to_attribute(&ctx, point_highlight_attr, 1 as i32, scene.points.highlights.as_slice());
|
||||
bind_new_buffer_to_attribute(&ctx, point_selection_attr, 1 as i32, scene.points.selections.as_slice());
|
||||
bind_new_buffer_to_attribute(&ctx, point_position_attr,
|
||||
SPACE_DIM as i32, point_positions.as_slice());
|
||||
bind_new_buffer_to_attribute(&ctx, point_color_attr,
|
||||
(COLOR_SIZE + 1) as i32,
|
||||
scene.points.colors_with_opacity.concat().as_slice());
|
||||
bind_new_buffer_to_attribute(&ctx, point_highlight_attr,
|
||||
1i32, scene.points.highlights.as_slice());
|
||||
bind_new_buffer_to_attribute(&ctx, point_selection_attr,
|
||||
1i32, scene.points.selections.as_slice());
|
||||
|
||||
// draw the scene
|
||||
ctx.draw_arrays(WebGl2RenderingContext::POINTS, 0, point_positions.ncols() as i32);
|
||||
ctx.draw_arrays(WebGl2RenderingContext::POINTS, 0,
|
||||
point_positions.ncols() as i32);
|
||||
|
||||
// disable the point program's vertex attributes
|
||||
ctx.disable_vertex_attrib_array(point_position_attr);
|
||||
|
@ -915,7 +957,9 @@ pub fn Display() -> View {
|
|||
.into_iter()
|
||||
.filter(|elt| !elt.ghost().get());
|
||||
for elt in tangible_elts {
|
||||
match assembly_to_world.with(|asm_to_world| elt.cast(dir, asm_to_world, pixel_size)) {
|
||||
let target = assembly_to_world.with(
|
||||
|asm_to_world| elt.cast(dir, asm_to_world, pixel_size));
|
||||
match target {
|
||||
Some(depth) => match clicked {
|
||||
Some((_, best_depth)) => {
|
||||
if depth < best_depth {
|
||||
|
|
|
@ -9,6 +9,7 @@ use crate::{
|
|||
Element,
|
||||
HalfCurvatureRegulator,
|
||||
InversiveDistanceRegulator,
|
||||
PointCoordinateRegulator,
|
||||
Regulator,
|
||||
},
|
||||
specified::SpecifiedValue
|
||||
|
@ -62,8 +63,10 @@ fn RegulatorInput(regulator: Rc<dyn Regulator>) -> View {
|
|||
placeholder = measurement.with(|result| result.to_string()),
|
||||
bind:value = value,
|
||||
on:change = move |_| {
|
||||
let specification =
|
||||
SpecifiedValue::try_from(value.get_clone_untracked());
|
||||
valid.set(
|
||||
match SpecifiedValue::try_from(value.get_clone_untracked()) {
|
||||
match specification {
|
||||
Ok(set_pt) => {
|
||||
set_point.set(set_pt);
|
||||
true
|
||||
|
@ -119,6 +122,20 @@ impl OutlineItem for HalfCurvatureRegulator {
|
|||
}
|
||||
}
|
||||
|
||||
impl OutlineItem for PointCoordinateRegulator {
|
||||
fn outline_item(self: Rc<Self>, _element: &Rc<dyn Element>) -> View {
|
||||
let name = format!("{} coordinate", self.axis);
|
||||
view! {
|
||||
li(class = "regulator") {
|
||||
div(class = "regulator-label") // for spacing
|
||||
div(class = "regulator-type") { (name) }
|
||||
RegulatorInput(regulator = self)
|
||||
div(class = "status")
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// a list item that shows an element in an outline view of an assembly
|
||||
#[component(inline_props)]
|
||||
fn ElementOutlineItem(element: Rc<dyn Element>) -> View {
|
||||
|
@ -126,7 +143,9 @@ fn ElementOutlineItem(element: Rc<dyn Element>) -> View {
|
|||
let class = {
|
||||
let element_for_class = element.clone();
|
||||
state.selection.map(
|
||||
move |sel| if sel.contains(&element_for_class) { "selected" } else { "" }
|
||||
move |sel|
|
||||
if sel.contains(&element_for_class) { "selected" }
|
||||
else { "" }
|
||||
)
|
||||
};
|
||||
let label = element.label().clone();
|
||||
|
@ -160,7 +179,8 @@ fn ElementOutlineItem(element: Rc<dyn Element>) -> View {
|
|||
move |event: KeyboardEvent| {
|
||||
match event.key().as_str() {
|
||||
"Enter" => {
|
||||
state.select(&element_for_handler, event.shift_key());
|
||||
state.select(&element_for_handler,
|
||||
event.shift_key());
|
||||
event.prevent_default();
|
||||
},
|
||||
"ArrowRight" if regulated.get() => {
|
||||
|
@ -190,18 +210,22 @@ fn ElementOutlineItem(element: Rc<dyn Element>) -> View {
|
|||
let state_for_handler = state.clone();
|
||||
let element_for_handler = element.clone();
|
||||
move |event: MouseEvent| {
|
||||
state_for_handler.select(&element_for_handler, event.shift_key());
|
||||
state_for_handler.select(&element_for_handler,
|
||||
event.shift_key());
|
||||
event.stop_propagation();
|
||||
event.prevent_default();
|
||||
}
|
||||
}
|
||||
) {
|
||||
div(class = "element-label") { (label) }
|
||||
div(class = "element-representation") { (rep_components) }
|
||||
div(class = "element-representation") {
|
||||
(rep_components)
|
||||
}
|
||||
input(
|
||||
r#type = "checkbox",
|
||||
bind:checked = element.ghost(),
|
||||
on:click = |event: MouseEvent| event.stop_propagation()
|
||||
on:click =
|
||||
|event: MouseEvent| event.stop_propagation()
|
||||
)
|
||||
}
|
||||
}
|
||||
|
|
|
@ -175,8 +175,9 @@ void main() {
|
|||
if (debug_mode) {
|
||||
// at the bottom of the screen, show the color scale instead of the
|
||||
// layer count
|
||||
if (gl_FragCoord.y < 10.) layer_cnt = int(16. * gl_FragCoord.x / resolution.x);
|
||||
|
||||
if (gl_FragCoord.y < 10.) {
|
||||
layer_cnt = int(16. * gl_FragCoord.x / resolution.x);
|
||||
}
|
||||
// convert number to color
|
||||
ivec3 bits = layer_cnt / ivec3(1, 2, 4);
|
||||
vec3 color = mod(vec3(bits), 2.);
|
||||
|
@ -217,14 +218,17 @@ void main() {
|
|||
// highlight intersections
|
||||
float ixn_dist = intersection_dist(frag, frag_next);
|
||||
float max_highlight = max(highlight, highlight_next);
|
||||
float ixn_highlight = 0.5 * max_highlight * (1. - smoothstep(2./3.*ixn_threshold, 1.5*ixn_threshold, ixn_dist));
|
||||
float ixn_highlight = 0.5 * max_highlight * (1. - smoothstep(
|
||||
2./3.*ixn_threshold, 1.5*ixn_threshold, ixn_dist));
|
||||
frag.color = mix(frag.color, vec4(1.), ixn_highlight);
|
||||
frag_next.color = mix(frag_next.color, vec4(1.), ixn_highlight);
|
||||
|
||||
// highlight cusps
|
||||
float cusp_cos = abs(dot(dir, frag.normal));
|
||||
float cusp_threshold = 2.*sqrt(ixn_threshold * sphere_list[hit.id].lt.s);
|
||||
float cusp_highlight = highlight * (1. - smoothstep(2./3.*cusp_threshold, 1.5*cusp_threshold, cusp_cos));
|
||||
float cusp_threshold = 2.*sqrt(
|
||||
ixn_threshold * sphere_list[hit.id].lt.s);
|
||||
float cusp_highlight = highlight * (1. - smoothstep(
|
||||
2./3.*cusp_threshold, 1.5*cusp_threshold, cusp_cos));
|
||||
frag.color = mix(frag.color, vec4(1.), cusp_highlight);
|
||||
|
||||
// composite the current fragment
|
||||
|
|
|
@ -167,29 +167,36 @@ fn load_low_curvature(assembly: &Assembly) {
|
|||
let curvature = plane.regulators().with_untracked(
|
||||
|regs| regs.first().unwrap().clone()
|
||||
);
|
||||
curvature.set_point().set(SpecifiedValue::try_from("0".to_string()).unwrap());
|
||||
curvature.set_point().set(
|
||||
SpecifiedValue::try_from("0".to_string()).unwrap());
|
||||
}
|
||||
let all_perpendicular = [central.clone()].into_iter()
|
||||
.chain(sides.clone())
|
||||
.chain(corners.clone());
|
||||
for sphere in all_perpendicular {
|
||||
// make each side and packed sphere perpendicular to the assembly plane
|
||||
let right_angle = InversiveDistanceRegulator::new([sphere, assemb_plane.clone()]);
|
||||
right_angle.set_point.set(SpecifiedValue::try_from("0".to_string()).unwrap());
|
||||
let right_angle = InversiveDistanceRegulator::new(
|
||||
[sphere, assemb_plane.clone()]);
|
||||
right_angle.set_point.set(
|
||||
SpecifiedValue::try_from("0".to_string()).unwrap());
|
||||
assembly.insert_regulator(Rc::new(right_angle));
|
||||
}
|
||||
for sphere in sides.clone().chain(corners.clone()) {
|
||||
// make each side and corner sphere tangent to the central sphere
|
||||
let tangency = InversiveDistanceRegulator::new([sphere.clone(), central.clone()]);
|
||||
tangency.set_point.set(SpecifiedValue::try_from("-1".to_string()).unwrap());
|
||||
let tangency = InversiveDistanceRegulator::new(
|
||||
[sphere.clone(), central.clone()]);
|
||||
tangency.set_point.set(
|
||||
SpecifiedValue::try_from("-1".to_string()).unwrap());
|
||||
assembly.insert_regulator(Rc::new(tangency));
|
||||
}
|
||||
for (side_index, side) in sides.enumerate() {
|
||||
// make each side tangent to the two adjacent corner spheres
|
||||
for (corner_index, corner) in corners.clone().enumerate() {
|
||||
if side_index != corner_index {
|
||||
let tangency = InversiveDistanceRegulator::new([side.clone(), corner]);
|
||||
tangency.set_point.set(SpecifiedValue::try_from("-1".to_string()).unwrap());
|
||||
let tangency = InversiveDistanceRegulator::new(
|
||||
[side.clone(), corner]);
|
||||
tangency.set_point.set(
|
||||
SpecifiedValue::try_from("-1".to_string()).unwrap());
|
||||
assembly.insert_regulator(Rc::new(tangency));
|
||||
}
|
||||
}
|
||||
|
@ -217,12 +224,15 @@ fn load_pointed(assembly: &Assembly) {
|
|||
for index_y in 0..=1 {
|
||||
let x = index_x as f64 - 0.5;
|
||||
let y = index_y as f64 - 0.5;
|
||||
|
||||
let x32 = x as f32;
|
||||
let y32 = y as f32;
|
||||
let coords =
|
||||
[0.5*(1.0 + x32), 0.5*(1.0 + y32), 0.5*(1.0 - x32*y32)];
|
||||
let _ = assembly.try_insert_element(
|
||||
Sphere::new(
|
||||
format!("sphere{index_x}{index_y}"),
|
||||
format!("Sphere {index_x}{index_y}"),
|
||||
[0.5*(1.0 + x) as f32, 0.5*(1.0 + y) as f32, 0.5*(1.0 - x*y) as f32],
|
||||
coords,
|
||||
engine::sphere(x, y, 0.0, 1.0),
|
||||
)
|
||||
);
|
||||
|
@ -231,7 +241,7 @@ fn load_pointed(assembly: &Assembly) {
|
|||
Point::new(
|
||||
format!("point{index_x}{index_y}"),
|
||||
format!("Point {index_x}{index_y}"),
|
||||
[0.5*(1.0 + x) as f32, 0.5*(1.0 + y) as f32, 0.5*(1.0 - x*y) as f32],
|
||||
coords,
|
||||
engine::point(x, y, 0.0),
|
||||
)
|
||||
);
|
||||
|
@ -320,19 +330,25 @@ fn load_tridiminished_icosahedron(assembly: &Assembly) {
|
|||
"face1".to_string(),
|
||||
"Face 1".to_string(),
|
||||
COLOR_FACE,
|
||||
engine::sphere_with_offset(frac_2_sqrt_6, -frac_1_sqrt_6, -frac_1_sqrt_6, -frac_1_sqrt_6, 0.0),
|
||||
engine::sphere_with_offset(
|
||||
frac_2_sqrt_6, -frac_1_sqrt_6, -frac_1_sqrt_6,
|
||||
-frac_1_sqrt_6, 0.0),
|
||||
),
|
||||
Sphere::new(
|
||||
"face2".to_string(),
|
||||
"Face 2".to_string(),
|
||||
COLOR_FACE,
|
||||
engine::sphere_with_offset(-frac_1_sqrt_6, frac_2_sqrt_6, -frac_1_sqrt_6, -frac_1_sqrt_6, 0.0),
|
||||
engine::sphere_with_offset(
|
||||
-frac_1_sqrt_6, frac_2_sqrt_6, -frac_1_sqrt_6,
|
||||
-frac_1_sqrt_6, 0.0),
|
||||
),
|
||||
Sphere::new(
|
||||
"face3".to_string(),
|
||||
"Face 3".to_string(),
|
||||
COLOR_FACE,
|
||||
engine::sphere_with_offset(-frac_1_sqrt_6, -frac_1_sqrt_6, frac_2_sqrt_6, -frac_1_sqrt_6, 0.0),
|
||||
engine::sphere_with_offset(
|
||||
-frac_1_sqrt_6, -frac_1_sqrt_6, frac_2_sqrt_6,
|
||||
-frac_1_sqrt_6, 0.0),
|
||||
),
|
||||
];
|
||||
for face in faces {
|
||||
|
@ -357,8 +373,10 @@ fn load_tridiminished_icosahedron(assembly: &Assembly) {
|
|||
let vertex_a = assembly.elements_by_id.with_untracked(
|
||||
|elts_by_id| elts_by_id[&format!("a{j}")].clone()
|
||||
);
|
||||
let incidence_a = InversiveDistanceRegulator::new([face.clone(), vertex_a.clone()]);
|
||||
incidence_a.set_point.set(SpecifiedValue::try_from("0".to_string()).unwrap());
|
||||
let incidence_a = InversiveDistanceRegulator::new(
|
||||
[face.clone(), vertex_a.clone()]);
|
||||
incidence_a.set_point.set(
|
||||
SpecifiedValue::try_from("0".to_string()).unwrap());
|
||||
assembly.insert_regulator(Rc::new(incidence_a));
|
||||
|
||||
// regulate the B-C vertex distances
|
||||
|
@ -380,13 +398,16 @@ fn load_tridiminished_icosahedron(assembly: &Assembly) {
|
|||
let vertex = assembly.elements_by_id.with_untracked(
|
||||
|elts_by_id| elts_by_id[&format!("{series}{k}")].clone()
|
||||
);
|
||||
let incidence = InversiveDistanceRegulator::new([face.clone(), vertex.clone()]);
|
||||
incidence.set_point.set(SpecifiedValue::try_from("0".to_string()).unwrap());
|
||||
let incidence = InversiveDistanceRegulator::new(
|
||||
[face.clone(), vertex.clone()]);
|
||||
incidence.set_point.set(
|
||||
SpecifiedValue::try_from("0".to_string()).unwrap());
|
||||
assembly.insert_regulator(Rc::new(incidence));
|
||||
|
||||
// regulate the A-B and A-C vertex distances
|
||||
assembly.insert_regulator(
|
||||
Rc::new(InversiveDistanceRegulator::new([vertex_a.clone(), vertex]))
|
||||
Rc::new(InversiveDistanceRegulator::new(
|
||||
[vertex_a.clone(), vertex]))
|
||||
);
|
||||
}
|
||||
}
|
||||
|
@ -434,7 +455,8 @@ fn load_dodecahedral_packing(assembly: &Assembly) {
|
|||
const COLOR_A: ElementColor = [1.00_f32, 0.25_f32, 0.00_f32];
|
||||
const COLOR_B: ElementColor = [1.00_f32, 0.00_f32, 0.25_f32];
|
||||
const COLOR_C: ElementColor = [0.25_f32, 0.00_f32, 1.00_f32];
|
||||
let phi = 0.5 + 1.25_f64.sqrt(); /* TO DO */ // replace with std::f64::consts::PHI when that gets stabilized
|
||||
/* TO DO */ // replace with std::f64::consts::PHI when that gets stabilized
|
||||
let phi = 0.5 + 1.25_f64.sqrt();
|
||||
let phi_inv = 1.0 / phi;
|
||||
let coord_scale = (phi + 2.0).sqrt();
|
||||
let face_scales = [phi_inv, (13.0 / 12.0) / coord_scale];
|
||||
|
@ -501,13 +523,16 @@ fn load_dodecahedral_packing(assembly: &Assembly) {
|
|||
|
||||
// make each face sphere perpendicular to the substrate
|
||||
for face in faces {
|
||||
let right_angle = InversiveDistanceRegulator::new([face, substrate.clone()]);
|
||||
right_angle.set_point.set(SpecifiedValue::try_from("0".to_string()).unwrap());
|
||||
let right_angle = InversiveDistanceRegulator::new(
|
||||
[face, substrate.clone()]);
|
||||
right_angle.set_point.set(
|
||||
SpecifiedValue::try_from("0".to_string()).unwrap());
|
||||
assembly.insert_regulator(Rc::new(right_angle));
|
||||
}
|
||||
|
||||
// set up the tangencies that define the packing
|
||||
for [long_edge_plane, short_edge_plane] in [["a", "b"], ["b", "c"], ["c", "a"]] {
|
||||
for [long_edge_plane, short_edge_plane]
|
||||
in [["a", "b"], ["b", "c"], ["c", "a"]] {
|
||||
for k in 0..2 {
|
||||
let long_edge_ids = [
|
||||
format!("{long_edge_plane}{k}0"),
|
||||
|
@ -526,9 +551,11 @@ fn load_dodecahedral_packing(assembly: &Assembly) {
|
|||
);
|
||||
|
||||
// set up the short-edge tangency
|
||||
let short_tangency = InversiveDistanceRegulator::new(short_edge.clone());
|
||||
let short_tangency = InversiveDistanceRegulator::new(
|
||||
short_edge.clone());
|
||||
if k == 0 {
|
||||
short_tangency.set_point.set(SpecifiedValue::try_from("-1".to_string()).unwrap());
|
||||
short_tangency.set_point.set(
|
||||
SpecifiedValue::try_from("-1".to_string()).unwrap());
|
||||
}
|
||||
assembly.insert_regulator(Rc::new(short_tangency));
|
||||
|
||||
|
@ -539,7 +566,9 @@ fn load_dodecahedral_packing(assembly: &Assembly) {
|
|||
[long_edge[i].clone(), short_edge[j].clone()]
|
||||
);
|
||||
if i == 0 && k == 0 {
|
||||
side_tangency.set_point.set(SpecifiedValue::try_from("-1".to_string()).unwrap());
|
||||
side_tangency.set_point.set(
|
||||
SpecifiedValue::try_from("-1".to_string()).unwrap()
|
||||
);
|
||||
}
|
||||
assembly.insert_regulator(Rc::new(side_tangency));
|
||||
}
|
||||
|
@ -604,7 +633,8 @@ fn load_balanced(assembly: &Assembly) {
|
|||
// initial configuration deliberately violates these constraints
|
||||
for inner in [a, b] {
|
||||
let tangency = InversiveDistanceRegulator::new([outer.clone(), inner]);
|
||||
tangency.set_point.set(SpecifiedValue::try_from("1".to_string()).unwrap());
|
||||
tangency.set_point.set(
|
||||
SpecifiedValue::try_from("1".to_string()).unwrap());
|
||||
assembly.insert_regulator(Rc::new(tangency));
|
||||
}
|
||||
}
|
||||
|
@ -712,10 +742,14 @@ fn load_radius_ratio(assembly: &Assembly) {
|
|||
[0.25_f32, 0.00_f32, 1.00_f32],
|
||||
].into_iter(),
|
||||
[
|
||||
engine::sphere_with_offset(base_dir[0], base_dir[1], base_dir[2], offset, 0.0),
|
||||
engine::sphere_with_offset(base_dir[0], -base_dir[1], -base_dir[2], offset, 0.0),
|
||||
engine::sphere_with_offset(-base_dir[0], base_dir[1], -base_dir[2], offset, 0.0),
|
||||
engine::sphere_with_offset(-base_dir[0], -base_dir[1], base_dir[2], offset, 0.0),
|
||||
engine::sphere_with_offset(
|
||||
base_dir[0], base_dir[1], base_dir[2], offset, 0.0),
|
||||
engine::sphere_with_offset(
|
||||
base_dir[0], -base_dir[1], -base_dir[2], offset, 0.0),
|
||||
engine::sphere_with_offset(
|
||||
-base_dir[0], base_dir[1], -base_dir[2], offset, 0.0),
|
||||
engine::sphere_with_offset(
|
||||
-base_dir[0], -base_dir[1], base_dir[2], offset, 0.0),
|
||||
].into_iter()
|
||||
).map(
|
||||
|(k, color, representation)| {
|
||||
|
@ -765,8 +799,10 @@ fn load_radius_ratio(assembly: &Assembly) {
|
|||
}
|
||||
|
||||
// put the vertices on the faces
|
||||
let incidence_regulator = InversiveDistanceRegulator::new([face_j.clone(), vertex_k.clone()]);
|
||||
incidence_regulator.set_point.set(SpecifiedValue::try_from("0".to_string()).unwrap());
|
||||
let incidence_regulator = InversiveDistanceRegulator::new(
|
||||
[face_j.clone(), vertex_k.clone()]);
|
||||
incidence_regulator.set_point.set(
|
||||
SpecifiedValue::try_from("0".to_string()).unwrap());
|
||||
assembly.insert_regulator(Rc::new(incidence_regulator));
|
||||
}
|
||||
}
|
||||
|
@ -860,25 +896,33 @@ fn load_irisawa_hexlet(assembly: &Assembly) {
|
|||
|elts_by_id| elts_by_id[&format!("chain{k}")].clone()
|
||||
)
|
||||
);
|
||||
for (chain_sphere, chain_sphere_next) in chain.clone().zip(chain.cycle().skip(1)) {
|
||||
for (chain_sphere, chain_sphere_next)
|
||||
in chain.clone().zip(chain.cycle().skip(1)) {
|
||||
for (other_sphere, inversive_distance) in [
|
||||
(outer.clone(), "1"),
|
||||
(sun.clone(), "-1"),
|
||||
(moon.clone(), "-1"),
|
||||
(chain_sphere_next.clone(), "-1"),
|
||||
] {
|
||||
let tangency = InversiveDistanceRegulator::new([chain_sphere.clone(), other_sphere]);
|
||||
tangency.set_point.set(SpecifiedValue::try_from(inversive_distance.to_string()).unwrap());
|
||||
let tangency = InversiveDistanceRegulator::new(
|
||||
[chain_sphere.clone(), other_sphere]);
|
||||
tangency.set_point.set(
|
||||
SpecifiedValue::try_from(
|
||||
inversive_distance.to_string()).unwrap());
|
||||
assembly.insert_regulator(Rc::new(tangency));
|
||||
}
|
||||
}
|
||||
|
||||
let outer_sun_tangency = InversiveDistanceRegulator::new([outer.clone(), sun]);
|
||||
outer_sun_tangency.set_point.set(SpecifiedValue::try_from("1".to_string()).unwrap());
|
||||
let outer_sun_tangency = InversiveDistanceRegulator::new(
|
||||
[outer.clone(), sun]);
|
||||
outer_sun_tangency.set_point.set(
|
||||
SpecifiedValue::try_from("1".to_string()).unwrap());
|
||||
assembly.insert_regulator(Rc::new(outer_sun_tangency));
|
||||
|
||||
let outer_moon_tangency = InversiveDistanceRegulator::new([outer.clone(), moon]);
|
||||
outer_moon_tangency.set_point.set(SpecifiedValue::try_from("1".to_string()).unwrap());
|
||||
let outer_moon_tangency = InversiveDistanceRegulator::new(
|
||||
[outer.clone(), moon]);
|
||||
outer_moon_tangency.set_point.set(
|
||||
SpecifiedValue::try_from("1".to_string()).unwrap());
|
||||
assembly.insert_regulator(Rc::new(outer_moon_tangency));
|
||||
}
|
||||
|
||||
|
@ -912,7 +956,8 @@ pub fn TestAssemblyChooser() -> View {
|
|||
"general" => load_general(assembly),
|
||||
"low-curvature" => load_low_curvature(assembly),
|
||||
"pointed" => load_pointed(assembly),
|
||||
"tridiminished-icosahedron" => load_tridiminished_icosahedron(assembly),
|
||||
"tridiminished-icosahedron" =>
|
||||
load_tridiminished_icosahedron(assembly),
|
||||
"dodecahedral-packing" => load_dodecahedral_packing(assembly),
|
||||
"balanced" => load_balanced(assembly),
|
||||
"off-center" => load_off_center(assembly),
|
||||
|
@ -929,7 +974,9 @@ pub fn TestAssemblyChooser() -> View {
|
|||
option(value = "general") { "General" }
|
||||
option(value = "low-curvature") { "Low-curvature" }
|
||||
option(value = "pointed") { "Pointed" }
|
||||
option(value = "tridiminished-icosahedron") { "Tridiminished icosahedron" }
|
||||
option(value = "tridiminished-icosahedron") {
|
||||
"Tridiminished icosahedron"
|
||||
}
|
||||
option(value = "dodecahedral-packing") { "Dodecahedral packing" }
|
||||
option(value = "balanced") { "Balanced" }
|
||||
option(value = "off-center") { "Off-center" }
|
||||
|
|
|
@ -9,8 +9,11 @@ pub fn point(x: f64, y: f64, z: f64) -> DVector<f64> {
|
|||
}
|
||||
|
||||
// the sphere with the given center and radius, with inward-pointing normals
|
||||
pub fn sphere(center_x: f64, center_y: f64, center_z: f64, radius: f64) -> DVector<f64> {
|
||||
let center_norm_sq = center_x * center_x + center_y * center_y + center_z * center_z;
|
||||
pub fn sphere(center_x: f64, center_y: f64, center_z: f64, radius: f64)
|
||||
-> DVector<f64>
|
||||
{
|
||||
let center_norm_sq =
|
||||
center_x * center_x + center_y * center_y + center_z * center_z;
|
||||
DVector::from_column_slice(&[
|
||||
center_x / radius,
|
||||
center_y / radius,
|
||||
|
@ -23,7 +26,9 @@ pub fn sphere(center_x: f64, center_y: f64, center_z: f64, radius: f64) -> DVect
|
|||
// the sphere of curvature `curv` whose closest point to the origin has position
|
||||
// `off * dir` and normal `dir`, where `dir` is a unit vector. setting the
|
||||
// curvature to zero gives a plane
|
||||
pub fn sphere_with_offset(dir_x: f64, dir_y: f64, dir_z: f64, off: f64, curv: f64) -> DVector<f64> {
|
||||
pub fn sphere_with_offset(
|
||||
dir_x: f64, dir_y: f64, dir_z: f64, off: f64, curv: f64) -> DVector<f64>
|
||||
{
|
||||
let norm_sp = 1.0 + off * curv;
|
||||
DVector::from_column_slice(&[
|
||||
norm_sp * dir_x,
|
||||
|
@ -52,8 +57,8 @@ pub fn project_point_to_normalized(rep: &mut DVector<f64>) {
|
|||
// --- partial matrices ---
|
||||
|
||||
pub struct MatrixEntry {
|
||||
index: (usize, usize),
|
||||
value: f64,
|
||||
pub index: (usize, usize),
|
||||
pub value: f64,
|
||||
}
|
||||
|
||||
pub struct PartialMatrix(Vec<MatrixEntry>);
|
||||
|
@ -200,7 +205,9 @@ impl ConfigSubspace {
|
|||
// 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<f64>, column_index: usize) -> DMatrix<f64> {
|
||||
pub fn proj(&self, v: &DVectorView<f64>, column_index: usize)
|
||||
-> DMatrix<f64>
|
||||
{
|
||||
if self.dim() == 0 {
|
||||
const ELEMENT_DIM: usize = 5;
|
||||
DMatrix::zeros(ELEMENT_DIM, self.assembly_dim)
|
||||
|
@ -291,7 +298,9 @@ impl SearchState {
|
|||
}
|
||||
}
|
||||
|
||||
fn basis_matrix(index: (usize, usize), nrows: usize, ncols: usize) -> DMatrix<f64> {
|
||||
fn basis_matrix(index: (usize, usize), nrows: usize, ncols: usize)
|
||||
-> DMatrix<f64>
|
||||
{
|
||||
let mut result = DMatrix::<f64>::zeros(nrows, ncols);
|
||||
result[index] = 1.0;
|
||||
result
|
||||
|
@ -414,7 +423,8 @@ pub fn realize_gram(
|
|||
for _ in 0..max_descent_steps {
|
||||
// find the negative gradient of the loss function
|
||||
let neg_grad = 4.0 * &*Q * &state.config * &state.err_proj;
|
||||
let mut neg_grad_stacked = neg_grad.clone().reshape_generic(Dyn(total_dim), Const::<1>);
|
||||
let mut neg_grad_stacked =
|
||||
neg_grad.clone().reshape_generic(Dyn(total_dim), Const::<1>);
|
||||
history.neg_grad.push(neg_grad.clone());
|
||||
|
||||
// find the negative Hessian of the loss function
|
||||
|
@ -431,7 +441,8 @@ pub fn realize_gram(
|
|||
-&basis_mat * &state.err_proj
|
||||
+ &state.config * &neg_d_err_proj
|
||||
);
|
||||
hess_cols.push(deriv_grad.reshape_generic(Dyn(total_dim), Const::<1>));
|
||||
hess_cols.push(
|
||||
deriv_grad.reshape_generic(Dyn(total_dim), Const::<1>));
|
||||
}
|
||||
}
|
||||
hess = DMatrix::from_columns(hess_cols.as_slice());
|
||||
|
@ -440,7 +451,8 @@ pub fn realize_gram(
|
|||
let hess_eigvals = hess.symmetric_eigenvalues();
|
||||
let min_eigval = hess_eigvals.min();
|
||||
if min_eigval <= 0.0 {
|
||||
hess -= reg_scale * min_eigval * DMatrix::identity(total_dim, total_dim);
|
||||
hess -= reg_scale * min_eigval
|
||||
* DMatrix::identity(total_dim, total_dim);
|
||||
}
|
||||
history.hess_eigvals.push(hess_eigvals);
|
||||
|
||||
|
@ -477,7 +489,8 @@ pub fn realize_gram(
|
|||
},
|
||||
};
|
||||
let base_step_stacked = hess_cholesky.solve(&neg_grad_stacked);
|
||||
let base_step = base_step_stacked.reshape_generic(Dyn(element_dim), Dyn(assembly_dim));
|
||||
let base_step = base_step_stacked.reshape_generic(
|
||||
Dyn(element_dim), Dyn(assembly_dim));
|
||||
history.base_step.push(base_step.clone());
|
||||
|
||||
// use backtracking line search to find a better configuration
|
||||
|
@ -507,9 +520,12 @@ pub fn realize_gram(
|
|||
}
|
||||
|
||||
// find the kernel of the Hessian. give it the uniform inner product
|
||||
let tangent = ConfigSubspace::symmetric_kernel(hess, unif_to_std, assembly_dim);
|
||||
let tangent =
|
||||
ConfigSubspace::symmetric_kernel(hess, unif_to_std, assembly_dim);
|
||||
|
||||
Ok(ConfigNeighborhood { #[cfg(feature = "dev")] config: state.config, nbhd: tangent })
|
||||
Ok(ConfigNeighborhood {
|
||||
#[cfg(feature = "dev")] config: state.config, nbhd: tangent
|
||||
})
|
||||
} else {
|
||||
Err("Failed to reach target accuracy".to_string())
|
||||
};
|
||||
|
@ -608,7 +624,8 @@ pub mod examples {
|
|||
for j in 0..2 {
|
||||
// diagonal and hinge edges
|
||||
for k in j..2 {
|
||||
problem.gram.push_sym(block + j, block + k, if j == k { 0.0 } else { -0.5 });
|
||||
problem.gram.push_sym(
|
||||
block + j, block + k, if j == k { 0.0 } else { -0.5 });
|
||||
}
|
||||
|
||||
// non-hinge edges
|
||||
|
@ -702,7 +719,8 @@ mod tests {
|
|||
]);
|
||||
for j in 0..2 {
|
||||
for k in j..2 {
|
||||
problem.gram.push_sym(j, k, if (j, k) == (1, 1) { 1.0 } else { 0.0 });
|
||||
problem.gram.push_sym(
|
||||
j, k, if (j, k) == (1, 1) { 1.0 } else { 0.0 });
|
||||
}
|
||||
}
|
||||
problem.frozen.push(3, 0, problem.guess[(3, 0)]);
|
||||
|
@ -729,7 +747,8 @@ mod tests {
|
|||
|
||||
// check against Irisawa's solution
|
||||
let entry_tol = SCALED_TOL.sqrt();
|
||||
let solution_diams = [30.0, 10.0, 6.0, 5.0, 15.0, 10.0, 3.75, 2.5, 2.0 + 8.0/11.0];
|
||||
let solution_diams =
|
||||
[30.0, 10.0, 6.0, 5.0, 15.0, 10.0, 3.75, 2.5, 2.0 + 8.0/11.0];
|
||||
for (k, diam) in solution_diams.into_iter().enumerate() {
|
||||
assert!((config[(3, k)] - 1.0 / diam).abs() < entry_tol);
|
||||
}
|
||||
|
@ -794,22 +813,29 @@ mod tests {
|
|||
// 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.column_iter().enumerate().map(
|
||||
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.column_iter().enumerate().map(
|
||||
|(k, v)| tangent.proj(&v, k)
|
||||
).sum();
|
||||
assert!((motion_std - motion_proj).norm_squared() < tol_sq);
|
||||
}
|
||||
}
|
||||
|
||||
fn translation_motion_unif(vel: &Vector3<f64>, assembly_dim: usize) -> Vec<DVector<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>> {
|
||||
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));
|
||||
|
@ -840,9 +866,12 @@ mod tests {
|
|||
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), config.column_iter().collect()),
|
||||
rotation_motion_unif(&Vector3::new(0.0, 1.0, 0.0), config.column_iter().collect()),
|
||||
rotation_motion_unif(&Vector3::new(0.0, 0.0, 1.0), config.column_iter().collect()),
|
||||
rotation_motion_unif(
|
||||
&Vector3::new(1.0, 0.0, 0.0), config.column_iter().collect()),
|
||||
rotation_motion_unif(
|
||||
&Vector3::new(0.0, 1.0, 0.0), config.column_iter().collect()),
|
||||
rotation_motion_unif(
|
||||
&Vector3::new(0.0, 0.0, 1.0), config.column_iter().collect()),
|
||||
|
||||
// the twist motion. more precisely: a motion that keeps the center
|
||||
// of mass stationary and preserves the distances between the
|
||||
|
@ -859,8 +888,10 @@ mod tests {
|
|||
[
|
||||
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]),
|
||||
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<_>>(),
|
||||
|
@ -880,9 +911,12 @@ mod tests {
|
|||
// 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(
|
||||
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);
|
||||
|
@ -913,10 +947,10 @@ mod tests {
|
|||
problem_orig.gram.push_sym(0, 0, 1.0);
|
||||
problem_orig.gram.push_sym(1, 1, 1.0);
|
||||
problem_orig.gram.push_sym(0, 1, 0.5);
|
||||
let Realization { result: result_orig, history: history_orig } = realize_gram(
|
||||
&problem_orig, SCALED_TOL, 0.5, 0.9, 1.1, 200, 110
|
||||
);
|
||||
let ConfigNeighborhood { config: config_orig, nbhd: tangent_orig } = result_orig.unwrap();
|
||||
let Realization { result: result_orig, history: history_orig } =
|
||||
realize_gram(&problem_orig, SCALED_TOL, 0.5, 0.9, 1.1, 200, 110);
|
||||
let ConfigNeighborhood { config: config_orig, nbhd: tangent_orig } =
|
||||
result_orig.unwrap();
|
||||
assert_eq!(config_orig, problem_orig.guess);
|
||||
assert_eq!(history_orig.scaled_loss.len(), 1);
|
||||
|
||||
|
@ -934,10 +968,10 @@ mod tests {
|
|||
frozen: problem_orig.frozen,
|
||||
guess: guess_tfm,
|
||||
};
|
||||
let Realization { result: result_tfm, history: history_tfm } = realize_gram(
|
||||
&problem_tfm, SCALED_TOL, 0.5, 0.9, 1.1, 200, 110
|
||||
);
|
||||
let ConfigNeighborhood { config: config_tfm, nbhd: tangent_tfm } = result_tfm.unwrap();
|
||||
let Realization { result: result_tfm, history: history_tfm } =
|
||||
realize_gram(&problem_tfm, SCALED_TOL, 0.5, 0.9, 1.1, 200, 110);
|
||||
let ConfigNeighborhood { config: config_tfm, nbhd: tangent_tfm } =
|
||||
result_tfm.unwrap();
|
||||
assert_eq!(config_tfm, problem_tfm.guess);
|
||||
assert_eq!(history_tfm.scaled_loss.len(), 1);
|
||||
|
||||
|
@ -948,7 +982,8 @@ mod tests {
|
|||
|
||||
// project the equivalent nudge to the tangent space of the solution
|
||||
// variety at the transformed solution
|
||||
let motion_tfm = DVector::from_column_slice(&[FRAC_1_SQRT_2, 0.0, FRAC_1_SQRT_2, 0.0]);
|
||||
let motion_tfm = DVector::from_column_slice(
|
||||
&[FRAC_1_SQRT_2, 0.0, FRAC_1_SQRT_2, 0.0]);
|
||||
let motion_tfm_proj = tangent_tfm.proj(&motion_tfm.as_view(), 0);
|
||||
|
||||
// take the transformation that sends the original solution to the
|
||||
|
@ -969,7 +1004,9 @@ mod tests {
|
|||
// the comparison tolerance because the transformation seems to
|
||||
// introduce some numerical error
|
||||
const SCALED_TOL_TFM: f64 = 1.0e-9;
|
||||
let tol_sq = ((problem_orig.guess.nrows() * problem_orig.guess.ncols()) as f64) * SCALED_TOL_TFM * SCALED_TOL_TFM;
|
||||
let tol_sq = ((problem_orig.guess.nrows()
|
||||
* problem_orig.guess.ncols()) as f64)
|
||||
* SCALED_TOL_TFM * SCALED_TOL_TFM;
|
||||
assert!((motion_proj_tfm - motion_tfm_proj).norm_squared() < tol_sq);
|
||||
}
|
||||
}
|
Loading…
Add table
Add a link
Reference in a new issue