Commit c93b8171 authored by SirBlueRabbit's avatar SirBlueRabbit
Browse files

maybe fixed 32 move

parent 3d774234
......@@ -17,8 +17,6 @@ fn main() -> Result<(), Box<dyn Error>> {
// precompute selection probabilities
let precomputed = Precomputed::new(&config.model);
// initialise step count
// equilibration phase
let mut step_count_therm = StepCount::default();
let mut triangulation = thermalise(&config, &precomputed, &mut step_count_therm);
......@@ -52,7 +50,7 @@ fn thermalise(
precomputed: &Precomputed,
step_count: &mut StepCount,
) -> Triangulation {
let mut triangulation = Triangulation::with_capacity(2 * config.model.sweep());
let mut triangulation = Triangulation::with_capacity(2 * config.model.sweep() + 12);
for _ in 0..(config.model.sweep() * config.markov_chain.thermalisation) {
triangulation.mc_step(config.model, precomputed, step_count);
}
......
......@@ -113,7 +113,7 @@ impl Triangulation {
dummy_label.next().next().adj_int(),
];
// replace old simplex with dummy
// glue dummy to old boundary
for i in 0..4 {
self.half_edges
.glue_triangles(dummy[i], boundary2[i].adj_ext(&self.half_edges));
......@@ -127,6 +127,9 @@ impl Triangulation {
// identify new boundary
let boundary3 = [b.adj_int(), c, a.adj_int(), c.adj_int()];
// update dummy adjacency
self.update_dummy_adj(dummy, boundary2, boundary3);
// replace dummy with new simplex
for i in 0..4 {
self.half_edges
......@@ -153,6 +156,80 @@ impl Triangulation {
fn move_32(&mut self, label: Label<HalfEdge>) {
// identify halfedges
let a = label
.next()
.next()
.adj_ext(&self.half_edges)
.adj_int()
.next()
.next();
let b = label
.adj_int()
.adj_ext(&self.half_edges)
.next()
.adj_int()
.next()
.next();
let c = label;
// identify vertices
let v0 = self.half_edges[c].vertex_tail;
let v1 = self.half_edges[a].vertex_tail;
let v2 = self.half_edges[a].vertex_head;
let v3 = self.half_edges[b].vertex_head;
let v4 = self.half_edges[c].vertex_head;
// create dummy simplex
let dummy_label = self.insert_tet([v0, v0, v0, v0]);
let dummy = [
dummy_label,
dummy_label.adj_int(),
dummy_label.next().adj_int(),
dummy_label.next().next().adj_int(),
];
// find old boundary
let bound3 = [
c.adj_int().next().next().adj_int().next(),
b,
c.next().next().adj_int().next(),
a,
];
// glue dummy to external boundary
for i in 0..4 {
self.half_edges
.glue_triangles(dummy[i], bound3[i].adj_ext(&self.half_edges));
}
// remove spare simplex
self.remove_tet(c);
// internal gluing
self.half_edges.glue_triangles(a, b.next().next());
// find triangles that should go to new boundary
let bound2 = [
a.next().adj_int(),
a.next().next().adj_int(),
b.next().adj_int(),
b.next().next().adj_int(),
];
// glue new boundary to dummy neighbours
for i in 0..4 {
self.half_edges
.glue_triangles(dummy[i].adj_ext(&self.half_edges), bound2[i]);
}
// remove dummy
self.remove_tet(dummy_label);
// update vertices
self.half_edges.set_vertices(a, [v1, v2, v3, v4]);
self.half_edges.set_vertices(b, [v1, v3, v2, v0]);
/* // identify halfedges
let a = label.next().next().adj_int().next();
let b = label
.adj_int()
......@@ -187,12 +264,16 @@ impl Triangulation {
dummy_label.next().next().adj_int(),
];
// replace old simplex with dummy
// glue dummy to old simplex neighbours
for i in 0..4 {
self.half_edges
.glue_triangles(dummy[i], boundary3[i].adj_ext(&self.half_edges));
}
// internal gluing
self.half_edges
.glue_triangles(a.adj_int(), b.adj_int().next().next());
// identify new boundary
let boundary2 = [
label.next().next().adj_int(),
......@@ -201,23 +282,38 @@ impl Triangulation {
label.adj_ext(&self.half_edges).next().next().adj_int(),
];
// replace dummy with new simplex
// update dummy adjacency
self.update_dummy_adj(dummy, boundary3, boundary2);
// glue dummy neighbours to new simplex
for i in 0..4 {
self.half_edges
.glue_triangles(dummy[i].adj_ext(&self.half_edges), boundary2[i]);
}
// remove extra tetrahedron and dummy
self.remove_tet(c);
// remove dummy and spare simplex
self.remove_tet(dummy_label);
self.remove_tet(c);
// internal gluing
self.half_edges
.glue_triangles(a.adj_int(), b.adj_int().next().next());
// match vertices
// update vertices
self.half_edges.set_vertices(a, [v1, v2, v0, v3]);
self.half_edges.set_vertices(b, [v3, v2, v4, v1]);
self.half_edges.set_vertices(b, [v3, v2, v4, v1]); */
}
fn update_dummy_adj(
&mut self,
dummy: [Label<HalfEdge>; 4],
old: [Label<HalfEdge>; 4],
new: [Label<HalfEdge>; 4],
) {
for i in 0..4 {
for j in 0..4 {
if dummy[i].adj_ext(&self.half_edges).same_triangle(old[j]) {
self.half_edges[dummy[i]].adj_ext =
Label::<HalfEdge>::from(new[j].value() / 3 * 3 + dummy[i].value() % 3);
}
}
}
}
/* fn shard_boundary(&self, label: Label<HalfEdge>, shard: Shard) -> Boundary {
......
......@@ -345,10 +345,10 @@ impl Label<HalfEdge> {
self.value() / 12 == other.value() / 12
}
/* pub fn same_triangle(&self, other: Label<HalfEdge>) -> bool {
pub fn same_triangle(&self, other: Label<HalfEdge>) -> bool {
self.value() / 3 == other.value() / 3
}
/*
pub fn rotation(&self, other: Label<HalfEdge>) -> Rotation {
match (self.value() + other.value()) % 3 {
0 => Rotation::Aligned,
......@@ -369,13 +369,25 @@ impl Pool<HalfEdge> {
self[b.next().next()].adj_ext = a.next();
}
pub fn set_vertices(&mut self, tet_label: Label<HalfEdge>, vertices: [Label<Vertex>; 4]) {
let index = tet_label.value() / 12 * 12;
(0..12).for_each(|i| {
let label = Label::<HalfEdge>::from(index + i);
self[label].vertex_tail = vertices[VERTICES[i][0]];
self[label].vertex_head = vertices[VERTICES[i][1]];
});
pub fn set_vertices(&mut self, label: Label<HalfEdge>, vertices: [Label<Vertex>; 4]) {
let tet_half_edges = [
label,
label.next(),
label.next().next(),
label.adj_int(),
label.adj_int().next(),
label.adj_int().next().next(),
label.next().adj_int().next().next(),
label.next().adj_int(),
label.next().adj_int().next(),
label.next().next().adj_int().next(),
label.next().next().adj_int().next().next(),
label.next().next().adj_int(),
];
for i in 0..12 {
self[tet_half_edges[i]].vertex_tail = vertices[VERTICES[i][0]];
self[tet_half_edges[i]].vertex_head = vertices[VERTICES[i][1]];
}
}
}
......@@ -482,7 +494,7 @@ mod tests {
// initialise
let model = Model::new(2, 0.01, 0.0, 3.0, 1, 1);
let precomputed = Precomputed::new(&model);
let mut triangulation = Triangulation::with_capacity(2 * model.sweep());
let mut triangulation = Triangulation::with_capacity(2 * model.sweep() + 12);
print!("{}", triangulation);
// perform a large number of steps and check each iteration
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment