forked from 1b-agent/fieldcraft-upstream
Compare commits
3
Commits
main
...
tet-84-speed
| Author | SHA1 | Date | |
|---|---|---|---|
 jyvikoandClaude Opus 5.5
|
e4290c7366 | ||
|
jyvikoandClaude Opus 5.5
|
4329630943 | ||
|
jyvikoandClaude Opus 5.5
|
22ef56f323 |
No files matched your search
@@ -275,7 +275,8 @@ fn on_box_edge(p: [f64; 3]) -> bool {
|
||||
fn strong_relaxation_of_a_filleted_box_keeps_faces_and_edges() {
|
||||
let (_d, mut s, _) = filleted_box();
|
||||
let before = s.model.clone();
|
||||
let r = s.execute("mesh.node.relax", &json!({"mode": "strong", "bodies": ["B"], "below": 0.3})).unwrap();
|
||||
// The tet mesher leaves no element below the skill's 0.3, so relax those below 0.45.
|
||||
let r = s.execute("mesh.node.relax", &json!({"mode": "strong", "bodies": ["B"], "below": 0.45})).unwrap();
|
||||
let (b, a) = (&r["before"], &r["after"]);
|
||||
// The worst tets may not improve: their corners can all be held on the box's edges.
|
||||
assert!(a["min_scaled_jacobian"].as_f64().unwrap() >= b["min_scaled_jacobian"].as_f64().unwrap(), "{r}");
|
||||
|
||||
@@ -170,7 +170,9 @@ pub enum Step {
|
||||
pub struct Raw {
|
||||
slots: Vec<(u32, [u32; 4], [u32; 4], [u32; 4])>,
|
||||
tets_len: usize,
|
||||
free: Vec<u32>,
|
||||
/// The free list's length, and its last entries (those the change may take).
|
||||
free_len: usize,
|
||||
free_tail: Vec<u32>,
|
||||
vtets: Vec<(u32, Vec<u32>)>,
|
||||
live: usize,
|
||||
}
|
||||
@@ -199,6 +201,21 @@ pub struct TetComplex {
|
||||
frozen: Vec<u32>,
|
||||
mark: Vec<u32>,
|
||||
epoch: u32,
|
||||
/// Change clock: `touched[v]` is the clock of the last change to a tetrahedron of `v`'s star
|
||||
/// (made, removed, reshaped or restored), `failed[t]` the clock at which every topological
|
||||
/// operation on slot `t` last failed ([`Self::known_to_fail`]).
|
||||
clock: u64,
|
||||
touched: Vec<u64>,
|
||||
failed: Vec<u64>,
|
||||
/// The clock at which removing an edge `(a, b)` (sorted), or a face from apex `a`'s side to apex `b` (sorted face)
|
||||
/// last failed: the operation reads only the tetrahedra holding both apexes.
|
||||
failed_edges: IdMap<(u32, u32), u64>,
|
||||
failed_faces: IdMap<(u32, u32, [u32; 3]), u64>,
|
||||
/// Per slot: the quality ([`super::measure::quality`]) and one more than the clock it was
|
||||
/// computed at; current while no vertex of the tetrahedron was touched since.
|
||||
quality_cache: std::cell::RefCell<Vec<(u64, f64)>>,
|
||||
/// The clock at which smoothing a vertex last found no move.
|
||||
failed_smoothing: Vec<u64>,
|
||||
live: usize,
|
||||
max_tets: usize,
|
||||
max_points: usize,
|
||||
@@ -281,6 +298,13 @@ impl TetComplex {
|
||||
frozen: Vec::new(),
|
||||
mark: vec![0; tets.len()],
|
||||
epoch: 0,
|
||||
clock: 0,
|
||||
touched: vec![0; points.len()],
|
||||
failed: Vec::new(),
|
||||
failed_edges: IdMap::default(),
|
||||
failed_faces: IdMap::default(),
|
||||
failed_smoothing: Vec::new(),
|
||||
quality_cache: std::cell::RefCell::new(Vec::new()),
|
||||
live: tets.len(),
|
||||
max_tets,
|
||||
max_points,
|
||||
@@ -355,6 +379,125 @@ impl TetComplex {
|
||||
if let Some(s) = self.kinds.get_mut(v as usize) {
|
||||
*s = k;
|
||||
}
|
||||
self.touch_star(v);
|
||||
}
|
||||
|
||||
// ---- change clock ----
|
||||
|
||||
fn touch(&mut self, vs: &[u32]) {
|
||||
self.clock += 1;
|
||||
for &v in vs {
|
||||
if let Some(s) = self.touched.get_mut(v as usize) {
|
||||
*s = self.clock;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/// Touch `v` and every vertex of its star (moving `v` reshapes those tetrahedra).
|
||||
fn touch_star(&mut self, v: u32) {
|
||||
let mut vs = vec![v];
|
||||
for &t in self.vtets.get(v as usize).map(Vec::as_slice).unwrap_or_default() {
|
||||
vs.extend(self.tets.get(t as usize).copied().unwrap_or([NONE; 4]));
|
||||
}
|
||||
self.touch(&vs);
|
||||
}
|
||||
|
||||
/// The last change near tetrahedron `t`: in the stars of its vertices and of its neighbours'
|
||||
/// vertices, the region every topological operation on `t` reads.
|
||||
fn stamp(&self, t: u32) -> Result<u64> {
|
||||
let mut m = 0;
|
||||
let mut look = |vs: [u32; 4]| {
|
||||
for v in vs {
|
||||
m = m.max(self.touched.get(v as usize).copied().unwrap_or(u64::MAX));
|
||||
}
|
||||
};
|
||||
look(self.tet(t)?);
|
||||
for n in self.adj(t)? {
|
||||
if n != NONE {
|
||||
look(self.tet(n)?);
|
||||
}
|
||||
}
|
||||
Ok(m)
|
||||
}
|
||||
|
||||
fn touched_at(&self, v: u32) -> u64 {
|
||||
self.touched.get(v as usize).copied().unwrap_or(u64::MAX)
|
||||
}
|
||||
|
||||
/// Record that removing edge `ab` failed now.
|
||||
pub fn mark_edge_failed(&mut self, a: u32, b: u32) {
|
||||
self.failed_edges.insert(sorted2(a, b), self.clock);
|
||||
}
|
||||
|
||||
/// Whether removing edge `ab` failed and no tetrahedron around `a` or `b` changed since.
|
||||
pub fn edge_known_to_fail(&self, a: u32, b: u32) -> bool {
|
||||
self.failed_edges.get(&sorted2(a, b)).is_some_and(|&f| self.touched_at(a).max(self.touched_at(b)) <= f)
|
||||
}
|
||||
|
||||
/// Record that removing `face` between apexes `a` and `b` failed now.
|
||||
pub fn mark_face_failed(&mut self, a: u32, b: u32, face: [u32; 3]) {
|
||||
self.failed_faces.insert((a, b, sorted3(face)), self.clock);
|
||||
}
|
||||
|
||||
/// Whether removing `face` between `a` and `b` failed and nothing around them changed since.
|
||||
pub fn face_known_to_fail(&self, a: u32, b: u32, face: [u32; 3]) -> bool {
|
||||
self.failed_faces.get(&(a, b, sorted3(face))).is_some_and(|&f| self.touched_at(a).max(self.touched_at(b)) <= f)
|
||||
}
|
||||
|
||||
/// The quality of live tetrahedron `t`, computed once per change of its vertices.
|
||||
pub fn slot_quality(&self, t: u32) -> Result<f64> {
|
||||
let vs = self.tet(t)?;
|
||||
let stamp = vs.iter().map(|&v| self.touched_at(v)).max().unwrap_or(u64::MAX);
|
||||
let i = t as usize;
|
||||
if let Some(&(c, q)) = self.quality_cache.borrow().get(i)
|
||||
&& c > stamp
|
||||
{
|
||||
return Ok(q);
|
||||
}
|
||||
let q = super::measure::quality(&self.coords(vs)?);
|
||||
let mut cache = self.quality_cache.borrow_mut();
|
||||
if cache.len() <= i {
|
||||
cache.resize(i + 1, (0, 0.0));
|
||||
}
|
||||
if let Some(e) = cache.get_mut(i) {
|
||||
*e = (self.clock + 1, q);
|
||||
}
|
||||
Ok(q)
|
||||
}
|
||||
|
||||
/// Record that smoothing `v` found no move now.
|
||||
pub fn mark_smoothing_failed(&mut self, v: u32) {
|
||||
let i = v as usize;
|
||||
if self.failed_smoothing.len() <= i {
|
||||
self.failed_smoothing.resize(i + 1, 0);
|
||||
}
|
||||
if let Some(f) = self.failed_smoothing.get_mut(i) {
|
||||
*f = self.clock;
|
||||
}
|
||||
}
|
||||
|
||||
/// Whether smoothing `v` found no move and no tetrahedron around it changed since.
|
||||
pub fn smoothing_known_to_fail(&self, v: u32) -> bool {
|
||||
self.failed_smoothing.get(v as usize).is_some_and(|&f| f > 0 && self.touched_at(v) <= f)
|
||||
}
|
||||
|
||||
/// Record that every topological operation on `t` failed now.
|
||||
pub fn mark_failed(&mut self, t: u32) {
|
||||
let i = t as usize;
|
||||
if self.failed.len() <= i {
|
||||
self.failed.resize(i + 1, 0);
|
||||
}
|
||||
if let Some(f) = self.failed.get_mut(i) {
|
||||
*f = self.clock;
|
||||
}
|
||||
}
|
||||
|
||||
/// Whether the topological operations on `t` failed and nothing they read has changed
|
||||
/// since, so trying again would fail the same way (they are deterministic functions of that
|
||||
/// region).
|
||||
pub fn known_to_fail(&self, t: u32) -> Result<bool> {
|
||||
let Some(&f) = self.failed.get(t as usize) else { return Ok(false) };
|
||||
Ok(f > 0 && self.stamp(t)? <= f)
|
||||
}
|
||||
|
||||
pub fn slots(&self) -> usize {
|
||||
@@ -424,6 +567,7 @@ impl TetComplex {
|
||||
|
||||
pub fn set_curve(&mut self, a: u32, b: u32, curve: Option<u32>) {
|
||||
let key = sorted2(a, b);
|
||||
self.touch(&[a, b]);
|
||||
let before = match curve {
|
||||
Some(c) => self.curves.insert(key, c),
|
||||
None => self.curves.remove(&key),
|
||||
@@ -688,6 +832,7 @@ impl TetComplex {
|
||||
self.pts.push(p);
|
||||
self.kinds.push(kind);
|
||||
self.vtets.push(Vec::new());
|
||||
self.touched.push(self.clock);
|
||||
if self.journal.recording {
|
||||
self.journal.steps.push(Step::Added { v });
|
||||
}
|
||||
@@ -695,7 +840,7 @@ impl TetComplex {
|
||||
}
|
||||
|
||||
/// The data a change of `slots` and of the incidence of `verts` would overwrite.
|
||||
fn snapshot(&self, slots: &[u32], verts: &[u32]) -> Raw {
|
||||
fn snapshot(&self, slots: &[u32], verts: &[u32], takes: usize) -> Raw {
|
||||
let mut ids: Vec<u32> = slots.to_vec();
|
||||
ids.sort_unstable();
|
||||
ids.dedup();
|
||||
@@ -711,7 +856,8 @@ impl TetComplex {
|
||||
})
|
||||
.collect(),
|
||||
tets_len: self.tets.len(),
|
||||
free: self.free.clone(),
|
||||
free_len: self.free.len(),
|
||||
free_tail: self.free.get(self.free.len().saturating_sub(takes)..).unwrap_or_default().to_vec(),
|
||||
vtets: vs.into_iter().filter_map(|v| Some((v, self.vtets.get(v as usize)?.clone()))).collect(),
|
||||
live: self.live,
|
||||
}
|
||||
@@ -719,6 +865,12 @@ impl TetComplex {
|
||||
|
||||
/// Put back what [`Self::snapshot`] saved.
|
||||
fn restore(&mut self, raw: Raw) {
|
||||
let mut touched: Vec<u32> = raw.vtets.iter().map(|(v, _)| *v).collect();
|
||||
for (t, vs, _, _) in &raw.slots {
|
||||
touched.extend(self.tets.get(*t as usize).copied().unwrap_or([NONE; 4]));
|
||||
touched.extend(vs);
|
||||
}
|
||||
self.touch(&touched);
|
||||
for (v, l) in raw.vtets {
|
||||
if let Some(s) = self.vtets.get_mut(v as usize) {
|
||||
*s = l;
|
||||
@@ -734,7 +886,8 @@ impl TetComplex {
|
||||
self.adj.truncate(raw.tets_len);
|
||||
self.tags.truncate(raw.tets_len);
|
||||
self.mark.truncate(raw.tets_len);
|
||||
self.free = raw.free;
|
||||
self.free.truncate(raw.free_len.saturating_sub(raw.free_tail.len()));
|
||||
self.free.extend(raw.free_tail);
|
||||
self.live = raw.live;
|
||||
}
|
||||
|
||||
@@ -744,6 +897,7 @@ impl TetComplex {
|
||||
self.pts.pop();
|
||||
self.kinds.pop();
|
||||
self.vtets.pop();
|
||||
self.touched.pop();
|
||||
return Ok(());
|
||||
}
|
||||
Err(corrupt("dropping a vertex in use"))
|
||||
@@ -755,6 +909,7 @@ impl TetComplex {
|
||||
let slot = self.pts.get_mut(v as usize).ok_or_else(|| corrupt("vertex id"))?;
|
||||
let from = *slot;
|
||||
*slot = p;
|
||||
self.touch_star(v);
|
||||
if self.journal.recording {
|
||||
self.journal.steps.push(Step::Moved { v, from });
|
||||
}
|
||||
@@ -962,7 +1117,7 @@ impl TetComplex {
|
||||
for &t in old {
|
||||
verts.extend(self.tet(t)?);
|
||||
}
|
||||
Some(Box::new(self.snapshot(&slots, &verts)))
|
||||
Some(Box::new(self.snapshot(&slots, &verts, change.new.len())))
|
||||
} else {
|
||||
None
|
||||
};
|
||||
@@ -1002,6 +1157,9 @@ impl TetComplex {
|
||||
}
|
||||
}
|
||||
}
|
||||
let mut changed: Vec<u32> = removed.iter().flat_map(|r| r.vs).collect();
|
||||
changed.extend(change.new.iter().flatten());
|
||||
self.touch(&changed);
|
||||
self.live = self.live - old.len() + ids.len();
|
||||
if self.journal.recording {
|
||||
self.journal.steps.push(Step::Replaced { removed, added: change.new.clone(), raw: raw.unwrap_or_default() });
|
||||
@@ -1044,7 +1202,7 @@ impl TetComplex {
|
||||
return Ok(false);
|
||||
}
|
||||
}
|
||||
let raw = self.journal.recording.then(|| Box::new(self.snapshot(&[t, ns[k], ns[l]], &vs)));
|
||||
let raw = self.journal.recording.then(|| Box::new(self.snapshot(&[t, ns[k], ns[l]], &vs, 0)));
|
||||
for m in [k, l] {
|
||||
let n = ns[m];
|
||||
let nn = self.adj(n)?;
|
||||
@@ -1055,6 +1213,7 @@ impl TetComplex {
|
||||
}
|
||||
}
|
||||
self.unlink(t, vs);
|
||||
self.touch(&vs);
|
||||
if let Some(s) = self.tets.get_mut(t as usize) {
|
||||
*s = [DEAD; 4];
|
||||
}
|
||||
@@ -1122,15 +1281,18 @@ impl TetComplex {
|
||||
if let Some(s) = self.pts.get_mut(v as usize) {
|
||||
*s = from;
|
||||
}
|
||||
self.touch_star(v);
|
||||
}
|
||||
Step::Added { v } => {
|
||||
if v as usize + 1 == self.pts.len() {
|
||||
self.pts.pop();
|
||||
self.kinds.pop();
|
||||
self.vtets.pop();
|
||||
self.touched.pop();
|
||||
}
|
||||
}
|
||||
Step::Curve { edge, curve } => {
|
||||
self.touch(&[edge.0, edge.1]);
|
||||
match curve {
|
||||
Some(c) => self.curves.insert(edge, c),
|
||||
None => self.curves.remove(&edge),
|
||||
|
||||
@@ -24,7 +24,7 @@ use super::complex::{EDGES, IdMap, NONE, Step, TetComplex, VertexKind, Walk, fac
|
||||
use super::measure::{self, P};
|
||||
use super::ops::{self, Site, remove_edge_at, remove_face, tet_quality};
|
||||
use super::refine::{Context, MAX_CAVITY, MIN_GAP};
|
||||
use super::smooth::smooth_vertex;
|
||||
use super::smooth::{smooth_vertex, smooth_vertex_below};
|
||||
|
||||
/// Settings of the schedule.
|
||||
#[derive(Debug, Clone, Copy, PartialEq)]
|
||||
@@ -53,6 +53,8 @@ pub struct ImproveStats {
|
||||
pub smoothed: usize,
|
||||
pub topological: usize,
|
||||
pub inserted: usize,
|
||||
/// Insertions tried (kept or undone).
|
||||
pub attempted: usize,
|
||||
/// The worst quality (smallest dihedral sine) after each pass.
|
||||
pub worst: Vec<f64>,
|
||||
pub stopped: Option<String>,
|
||||
@@ -81,8 +83,8 @@ fn indicators(cx: &TetComplex) -> Result<[f64; 8]> {
|
||||
let mut out = [0.0; 8];
|
||||
out[0] = f64::INFINITY;
|
||||
let mut n = 0usize;
|
||||
for (_, vs) in cx.live_tets() {
|
||||
let q = tet_quality(cx, vs)?;
|
||||
for (t, _) in cx.live_tets() {
|
||||
let q = cx.slot_quality(t)?;
|
||||
out[0] = out[0].min(q);
|
||||
for (i, s) in SIN_DEG.iter().enumerate() {
|
||||
out[i + 1] += q.min(*s);
|
||||
@@ -102,8 +104,8 @@ fn progressed(now: &[f64; 8], before: &[f64; 8]) -> bool {
|
||||
/// Smooth every movable vertex of a tetrahedron below `good`. Returns the moves.
|
||||
fn smoothing_pass(cx: &mut TetComplex, ctx: &Context, good: f64) -> Result<usize> {
|
||||
let mut want = vec![false; cx.point_count()];
|
||||
for (_, vs) in cx.live_tets() {
|
||||
if tet_quality(cx, vs)? < good {
|
||||
for (t, vs) in cx.live_tets() {
|
||||
if cx.slot_quality(t)? < good {
|
||||
for v in vs {
|
||||
if let Some(w) = want.get_mut(v as usize) {
|
||||
*w = true;
|
||||
@@ -149,7 +151,7 @@ fn delete_flat(cx: &mut TetComplex, t: u32) -> Result<bool> {
|
||||
/// edge ([`ops::collapse`]), which takes the tetrahedron away.
|
||||
fn collapse_into_crease(cx: &mut TetComplex, t: u32) -> Result<bool> {
|
||||
let vs = cx.tet(t)?;
|
||||
if tet_quality(cx, vs)? > 0.1 {
|
||||
if cx.slot_quality(t)? > 0.1 {
|
||||
return Ok(false);
|
||||
}
|
||||
let (ns, tags) = (cx.adj(t)?, cx.tags(t)?);
|
||||
@@ -182,7 +184,7 @@ fn collapse_into_crease(cx: &mut TetComplex, t: u32) -> Result<bool> {
|
||||
pub fn topological_pass(cx: &mut TetComplex, list: &[u32]) -> Result<usize> {
|
||||
let mut changes = 0;
|
||||
for &t in list {
|
||||
if !cx.is_live(t) {
|
||||
if !cx.is_live(t) || cx.known_to_fail(t)? {
|
||||
continue;
|
||||
}
|
||||
let vs = cx.tet(t)?;
|
||||
@@ -201,20 +203,25 @@ pub fn topological_pass(cx: &mut TetComplex, list: &[u32]) -> Result<usize> {
|
||||
changes += 1;
|
||||
continue;
|
||||
}
|
||||
let mut removed = false;
|
||||
for k in 0..4 {
|
||||
if remove_face(cx, t, k)? {
|
||||
changes += 1;
|
||||
removed = true;
|
||||
break;
|
||||
}
|
||||
}
|
||||
if !removed {
|
||||
cx.mark_failed(t);
|
||||
}
|
||||
}
|
||||
Ok(changes)
|
||||
}
|
||||
|
||||
fn below(cx: &TetComplex, good: f64) -> Result<Vec<u32>> {
|
||||
let mut out = Vec::new();
|
||||
for (t, vs) in cx.live_tets() {
|
||||
if tet_quality(cx, vs)? < good {
|
||||
for (t, _) in cx.live_tets() {
|
||||
if cx.slot_quality(t)? < good {
|
||||
out.push(t);
|
||||
}
|
||||
}
|
||||
@@ -288,7 +295,8 @@ fn change_quality(cx: &mut TetComplex, len: usize) -> Result<(f64, f64)> {
|
||||
}
|
||||
|
||||
/// The composite insertion of §3.5 (see the module docs). Returns whether it was kept.
|
||||
fn attempt_insert(cx: &mut TetComplex, ctx: &Context, p: Point3, kind: VertexKind, site: Site) -> Result<bool> {
|
||||
fn attempt_insert(cx: &mut TetComplex, ctx: &Context, attempted: &mut usize, good: f64, p: Point3, kind: VertexKind, site: Site) -> Result<bool> {
|
||||
*attempted += 1;
|
||||
let was = cx.journal.recording;
|
||||
cx.journal.recording = true;
|
||||
let len = cx.journal.steps.len();
|
||||
@@ -303,8 +311,15 @@ fn attempt_insert(cx: &mut TetComplex, ctx: &Context, p: Point3, kind: VertexKin
|
||||
if !cx.is_used(v) {
|
||||
break;
|
||||
}
|
||||
let star = cx.star(v)?;
|
||||
if topological_pass(cx, &star)? == 0 {
|
||||
// Only the tetrahedra below the threshold, as in the passes.
|
||||
let mut star = cx.star(v)?;
|
||||
let mut keep = Vec::with_capacity(star.len());
|
||||
for &t in &star {
|
||||
keep.push(cx.slot_quality(t)? < good);
|
||||
}
|
||||
let mut flags = keep.into_iter();
|
||||
star.retain(|_| flags.next().unwrap_or(false));
|
||||
if star.is_empty() || topological_pass(cx, &star)? == 0 {
|
||||
break;
|
||||
}
|
||||
}
|
||||
@@ -312,7 +327,7 @@ fn attempt_insert(cx: &mut TetComplex, ctx: &Context, p: Point3, kind: VertexKin
|
||||
let mut around: Vec<u32> = cx.neighbours(v)?;
|
||||
around.push(v);
|
||||
for u in around {
|
||||
smooth_vertex(cx, ctx, u)?;
|
||||
smooth_vertex_below(cx, ctx, u, good)?;
|
||||
}
|
||||
}
|
||||
let (before, after) = change_quality(cx, len)?;
|
||||
@@ -337,7 +352,7 @@ fn attempt_insert(cx: &mut TetComplex, ctx: &Context, p: Point3, kind: VertexKin
|
||||
}
|
||||
|
||||
/// Try inserting a vertex in each tetrahedron of `list` (§4.2). Returns the insertions kept.
|
||||
fn insertion_pass(cx: &mut TetComplex, ctx: &Context, list: &[u32], budget: usize, good: f64) -> Result<usize> {
|
||||
fn insertion_pass(cx: &mut TetComplex, ctx: &Context, list: &[u32], budget: usize, good: f64, attempted: &mut usize) -> Result<usize> {
|
||||
let mut kept = 0;
|
||||
for &t in list {
|
||||
if kept >= budget {
|
||||
@@ -347,7 +362,7 @@ fn insertion_pass(cx: &mut TetComplex, ctx: &Context, list: &[u32], budget: usiz
|
||||
continue;
|
||||
}
|
||||
let vs = cx.tet(t)?;
|
||||
if tet_quality(cx, vs)? >= good {
|
||||
if cx.slot_quality(t)? >= good {
|
||||
continue;
|
||||
}
|
||||
let ns = cx.adj(t)?;
|
||||
@@ -374,7 +389,7 @@ fn insertion_pass(cx: &mut TetComplex, ctx: &Context, list: &[u32], budget: usiz
|
||||
let target = measure::add(measure::mid(pa, pb), measure::scale(inward, 0.4 * l / len));
|
||||
if let Ok(p) = Point3::try_from(target)
|
||||
&& let Some(Walk::Inside(at)) = cx.walk(t, &p, 1 << 12)?
|
||||
&& attempt_insert(cx, ctx, p, VertexKind::Interior, Site::Interior(at))?
|
||||
&& attempt_insert(cx, ctx, attempted, good, p, VertexKind::Interior, Site::Interior(at))?
|
||||
{
|
||||
done = true;
|
||||
}
|
||||
@@ -391,7 +406,7 @@ fn insertion_pass(cx: &mut TetComplex, ctx: &Context, list: &[u32], budget: usiz
|
||||
let kind = VertexKind::Face(tags[k]);
|
||||
if let Some(p) = ctx.place(kind, g)
|
||||
&& cx.is_live(t)
|
||||
&& attempt_insert(cx, ctx, p, kind, Site::Face(t, k))?
|
||||
&& attempt_insert(cx, ctx, attempted, good, p, kind, Site::Face(t, k))?
|
||||
{
|
||||
done = true;
|
||||
}
|
||||
@@ -400,7 +415,7 @@ fn insertion_pass(cx: &mut TetComplex, ctx: &Context, list: &[u32], budget: usiz
|
||||
if !done && cx.is_live(t) && cx.tet(t)? == vs {
|
||||
let g = measure::centroid(&cx.coords(vs)?);
|
||||
if let Ok(p) = Point3::try_from(g)
|
||||
&& attempt_insert(cx, ctx, p, VertexKind::Interior, Site::Interior(t))?
|
||||
&& attempt_insert(cx, ctx, attempted, good, p, VertexKind::Interior, Site::Interior(t))?
|
||||
{
|
||||
done = true;
|
||||
}
|
||||
@@ -430,7 +445,7 @@ fn insertion_pass(cx: &mut TetComplex, ctx: &Context, list: &[u32], budget: usiz
|
||||
};
|
||||
let m = measure::mid(cx.xyz(a)?, cx.xyz(b)?);
|
||||
if let Some(p) = ctx.place(kind, m)
|
||||
&& attempt_insert(cx, ctx, p, kind, Site::Edge(a, b))?
|
||||
&& attempt_insert(cx, ctx, attempted, good, p, kind, Site::Edge(a, b))?
|
||||
{
|
||||
done = true;
|
||||
break;
|
||||
@@ -479,8 +494,8 @@ pub fn improve(cx: &mut TetComplex, ctx: &Context, o: &ImproveOptions) -> Result
|
||||
// The worst 3.5% of the tetrahedra, or (desperation pass, §4.2) every one with a
|
||||
// dihedral angle below 40° or above 140° (sine below sin 40°).
|
||||
let mut scored: Vec<(f64, u32)> = Vec::new();
|
||||
for (t, vs) in cx.live_tets() {
|
||||
scored.push((tet_quality(cx, vs)?, t));
|
||||
for (t, _) in cx.live_tets() {
|
||||
scored.push((cx.slot_quality(t)?, t));
|
||||
}
|
||||
scored.sort_by(|a, b| a.0.total_cmp(&b.0).then(a.1.cmp(&b.1)));
|
||||
let list: Vec<u32> = if failed == 1 {
|
||||
@@ -490,7 +505,7 @@ pub fn improve(cx: &mut TetComplex, ctx: &Context, o: &ImproveOptions) -> Result
|
||||
scored.iter().take(n).map(|s| s.1).collect()
|
||||
};
|
||||
let budget = o.max_insertions - stats.inserted;
|
||||
stats.inserted += insertion_pass(cx, ctx, &list, budget, o.good)?;
|
||||
stats.inserted += insertion_pass(cx, ctx, &list, budget, o.good, &mut stats.attempted)?;
|
||||
}
|
||||
let now = indicators(cx)?;
|
||||
carried = Some(now);
|
||||
@@ -548,7 +563,7 @@ mod tests {
|
||||
let q = Point3::try_from(q).unwrap();
|
||||
if let Some(Walk::Inside(at)) = cx.walk(t, &q, 1 << 12).unwrap() {
|
||||
tried += 1;
|
||||
assert!(!attempt_insert(&mut cx, &ctx, q, VertexKind::Interior, Site::Interior(at)).unwrap());
|
||||
assert!(!attempt_insert(&mut cx, &ctx, &mut 0, SIN_40, q, VertexKind::Interior, Site::Interior(at)).unwrap());
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
@@ -74,12 +74,41 @@ pub fn dihedral_cosines(x: &[P; 4]) -> [f64; 6] {
|
||||
})
|
||||
}
|
||||
|
||||
/// The quality objective: the smallest signed dihedral sine (1 is impossible; the regular
|
||||
/// tetrahedron scores √(8/9) ≈ 0.943). Positive for a valid tetrahedron, 0 for a flat one,
|
||||
/// negative for an inverted one (then the sine of largest magnitude, so worse inversions score
|
||||
/// lower). With `sin θ_ij = 6V ℓ_ij / (2A_k · 2A_l)`, the extreme sine is found on squared
|
||||
/// lengths and areas, with one square root.
|
||||
/// The quality objective of a valid tetrahedron: the smaller of its smallest dihedral sine and
|
||||
/// its weighted scaled Jacobian ([`JACOBIAN_WEIGHT`]; the regular tetrahedron scores
|
||||
/// √(8/9) ≈ 0.943). The dihedral sine
|
||||
/// catches slivers and caps; it does not see needles, whose dihedral angles can all be fine
|
||||
/// while a corner is nearly flat (Shewchuk 2002), and the scaled Jacobian
|
||||
/// `√2 · 6V / (ℓ₁ℓ₂ℓ₃)` at the worst corner, the measure analysts check (the skill's ≥ 0.2,
|
||||
/// [`crate::quality`]), does. 0 for a flat tetrahedron; for an inverted one the negative sine of
|
||||
/// largest magnitude, so worse inversions score lower. Extremes are found on squared lengths and
|
||||
/// areas, with one square root each.
|
||||
pub fn quality(x: &[P; 4]) -> f64 {
|
||||
quality_in_order(&canonical_order(x))
|
||||
}
|
||||
|
||||
/// The corners reordered by an even permutation that depends only on the corners' positions
|
||||
/// (sorted by coordinates, the last two swapped back when the sort was odd), so a measure
|
||||
/// computed on it gives the same bits for every listing of a tetrahedron with the same
|
||||
/// orientation: a removal tried from either side of a ring sees the same qualities.
|
||||
fn canonical_order(x: &[P; 4]) -> [P; 4] {
|
||||
let less = |a: &P, b: &P| a[0] < b[0] || (a[0] == b[0] && (a[1] < b[1] || (a[1] == b[1] && a[2] < b[2])));
|
||||
let mut y = *x;
|
||||
let mut odd = false;
|
||||
// A sorting network on four (pairs 01, 23, 02, 13, 12), counting swaps for the parity.
|
||||
for (i, j) in [(0, 1), (2, 3), (0, 2), (1, 3), (1, 2)] {
|
||||
if less(&y[j], &y[i]) {
|
||||
y.swap(i, j);
|
||||
odd = !odd;
|
||||
}
|
||||
}
|
||||
if odd {
|
||||
y.swap(2, 3);
|
||||
}
|
||||
y
|
||||
}
|
||||
|
||||
fn quality_in_order(x: &[P; 4]) -> f64 {
|
||||
let v6 = volume6(x);
|
||||
if v6 == 0.0 || !v6.is_finite() {
|
||||
return 0.0;
|
||||
@@ -95,9 +124,124 @@ pub fn quality(x: &[P; 4]) -> f64 {
|
||||
pick = if v6 > 0.0 { pick.min(r) } else { pick.max(r) };
|
||||
}
|
||||
let q = v6 * pick.sqrt();
|
||||
let q = if v6 > 0.0 { q.min(JACOBIAN_WEIGHT * scaled_jacobian(x, v6)) } else { q };
|
||||
if q.is_finite() { q } else { 0.0 }
|
||||
}
|
||||
|
||||
/// The corners of a tetrahedron and the other three corners of each.
|
||||
const CORNERS: [(usize, [usize; 3]); 4] = [(0, [1, 2, 3]), (1, [0, 2, 3]), (2, [0, 1, 3]), (3, [0, 1, 2])];
|
||||
|
||||
/// The squared product of the three edge lengths at each corner.
|
||||
fn corner_products(x: &[P; 4]) -> [f64; 4] {
|
||||
CORNERS.map(|(c, o)| {
|
||||
o.iter().fold(1.0, |m, &j| {
|
||||
let e = sub(x[j], x[c]);
|
||||
m * dot(e, e)
|
||||
})
|
||||
})
|
||||
}
|
||||
|
||||
/// The scaled Jacobian at the worst corner of a tetrahedron with six times the volume `v6`
|
||||
/// (as [`crate::quality`] reports it).
|
||||
fn scaled_jacobian(x: &[P; 4], v6: f64) -> f64 {
|
||||
let worst = corner_products(x).iter().fold(0.0_f64, |m, &p| m.max(p));
|
||||
std::f64::consts::SQRT_2 * v6 / worst.sqrt()
|
||||
}
|
||||
|
||||
/// The weight of the scaled Jacobian in [`quality`]: the meshing skill's limit, a scaled
|
||||
/// Jacobian of 0.3, counts as much as a 30° dihedral angle (sin 30°, the improvement's target,
|
||||
/// [`super::improve::ImproveOptions::target`]), so the Jacobian decides only for tetrahedra near
|
||||
/// that limit (needles) and the dihedral sine everywhere else.
|
||||
pub const JACOBIAN_WEIGHT: f64 = 0.5 / 0.3;
|
||||
|
||||
/// [`quality`] and its gradient with respect to corner `i`, for a positive tetrahedron: the
|
||||
/// gradient of the term [`quality`] picks (the smallest dihedral sine, or the scaled Jacobian at
|
||||
/// the worst corner, `∇ ln J = ∇V₆ / V₆ − Σ ∇ℓ / ℓ` over that corner's edges), analytic, as Klingner and
|
||||
/// Shewchuk (§3.1) and Freitag and Ollivier-Gooch use for smoothing. With
|
||||
/// `s = V₆ ℓ / (N_k N_l)` (`N` twice a face's area),
|
||||
/// `∇s = s (∇V₆ / V₆ + ∇ℓ / ℓ − ∇N_k / N_k − ∇N_l / N_l)`, where `∇V₆` at a corner is its
|
||||
/// opposite face's area vector, `∇ℓ` the unit edge vector at an end of the edge and, for a face
|
||||
/// `(p, q, r)` with unit normal `n̂`, `∇_p N = (q − r) × n̂`. A flat or inverted tetrahedron
|
||||
/// gets a zero gradient (smoothing only moves within valid positions).
|
||||
pub fn quality_gradient(x: &[P; 4], i: usize) -> (f64, P) {
|
||||
let q = quality(x);
|
||||
let v6 = volume6(x);
|
||||
if !(v6 > 0.0 && q > 0.0) || i > 3 {
|
||||
return (q, [0.0; 3]);
|
||||
}
|
||||
let crosses = [0, 1, 2, 3].map(|k| face_cross(x, k));
|
||||
let n2 = crosses.map(|c| dot(c, c));
|
||||
// The edge whose sine is smallest (the same comparison as `quality`).
|
||||
let mut pick = (f64::INFINITY, EDGE_FACES[0]);
|
||||
for e in EDGE_FACES {
|
||||
let (a, b, k, l) = e;
|
||||
let d = sub(x[a], x[b]);
|
||||
let r = dot(d, d) / (n2[k] * n2[l]);
|
||||
if r < pick.0 {
|
||||
pick = (r, e);
|
||||
}
|
||||
}
|
||||
// The scaled Jacobian branch, when it is the smaller.
|
||||
let products = corner_products(x);
|
||||
let (c, cp) = products.iter().enumerate().fold((0, 0.0_f64), |m, (c, &p)| if p > m.1 { (c, p) } else { m });
|
||||
let sj = JACOBIAN_WEIGHT * std::f64::consts::SQRT_2 * v6 / cp.sqrt();
|
||||
let sine = v6 * pick.0.sqrt();
|
||||
let face = |k: usize| -> [usize; 3] {
|
||||
match k {
|
||||
0 => [1, 2, 3],
|
||||
1 => [0, 2, 3],
|
||||
2 => [0, 1, 3],
|
||||
_ => [0, 1, 2],
|
||||
}
|
||||
};
|
||||
// ∇V₆ at corner i: the area vector of the face opposite i, signed so that V₆ grows along it.
|
||||
let gi = crosses[i];
|
||||
let gv = if dot(gi, sub(x[i], x[face(i)[0]])) < 0.0 { scale(gi, -1.0) } else { gi };
|
||||
if sj < sine {
|
||||
// ∇ ln sj = ∇V₆ / V₆ − Σ ∇ℓ / ℓ over the corner's edges at i.
|
||||
let (corner, others) = CORNERS[c];
|
||||
let mut g = scale(gv, 1.0 / v6);
|
||||
for j in others {
|
||||
let (p, o) = if i == corner {
|
||||
(x[i], x[j])
|
||||
} else if i == j {
|
||||
(x[i], x[corner])
|
||||
} else {
|
||||
continue;
|
||||
};
|
||||
let e = sub(p, o);
|
||||
g = sub(g, scale(e, 1.0 / dot(e, e)));
|
||||
}
|
||||
let g = scale(g, q);
|
||||
return if g.iter().all(|c| c.is_finite()) { (q, g) } else { (q, [0.0; 3]) };
|
||||
}
|
||||
let (a, b, k, l) = pick.1;
|
||||
let len = dist(x[a], x[b]);
|
||||
let gl = if i == a {
|
||||
scale(sub(x[a], x[b]), 1.0 / len)
|
||||
} else if i == b {
|
||||
scale(sub(x[b], x[a]), 1.0 / len)
|
||||
} else {
|
||||
[0.0; 3]
|
||||
};
|
||||
// ∇N of the face opposite m at corner i (zero when the face does not hold i).
|
||||
let gn = |m: usize| -> P {
|
||||
let f = face(m);
|
||||
let Some(pos) = f.iter().position(|&c| c == i) else { return [0.0; 3] };
|
||||
let n = crosses[m];
|
||||
let nn = norm(n);
|
||||
if !(nn > 0.0) {
|
||||
return [0.0; 3];
|
||||
}
|
||||
let (qv, rv) = (x[f[(pos + 1) % 3]], x[f[(pos + 2) % 3]]);
|
||||
scale(cross(sub(qv, rv), n), 1.0 / nn)
|
||||
};
|
||||
let (nk, nl) = (n2[k].sqrt(), n2[l].sqrt());
|
||||
let g = add(add(scale(gv, 1.0 / v6), scale(gl, 1.0 / len)), scale(add(scale(gn(k), 1.0 / nk), scale(gn(l), 1.0 / nl)), -1.0));
|
||||
let g = scale(g, q);
|
||||
if g.iter().all(|c| c.is_finite()) { (q, g) } else { (q, [0.0; 3]) }
|
||||
}
|
||||
|
||||
/// The smallest and largest dihedral angle in degrees (reports only).
|
||||
pub fn dihedral_range(x: &[P; 4]) -> (f64, f64) {
|
||||
let (s, c) = (dihedral_sines(x), dihedral_cosines(x));
|
||||
@@ -221,6 +365,49 @@ mod tests {
|
||||
assert!(circumcenter(&flat).is_none());
|
||||
}
|
||||
|
||||
/// Every even reordering of the corners gives the same quality, bit for bit.
|
||||
#[test]
|
||||
fn quality_does_not_depend_on_the_corner_order() {
|
||||
let x = positive([[0.1, 0.0, 0.05], [1.0, 0.2, 0.0], [0.3, 0.9, 0.1], [0.4, 0.35, 0.7]]);
|
||||
let q = quality(&x);
|
||||
for p in [[1, 2, 0, 3], [2, 0, 1, 3], [0, 2, 3, 1], [3, 1, 0, 2], [1, 0, 3, 2], [2, 3, 0, 1]] {
|
||||
assert_eq!(quality(&p.map(|i| x[i])).to_bits(), q.to_bits(), "{p:?}");
|
||||
}
|
||||
}
|
||||
|
||||
/// A needle: every dihedral angle above 20°, yet one corner nearly flat; the objective sees
|
||||
/// it through the scaled Jacobian.
|
||||
#[test]
|
||||
fn a_needle_scores_low() {
|
||||
let x = positive([[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.5, 0.866, 0.0], [0.5, 0.29, 12.0]]);
|
||||
let (lo, hi) = dihedral_range(&x);
|
||||
assert!(lo > 20.0 && hi < 120.0, "{lo} {hi}");
|
||||
assert!(quality(&x) < 0.15 * JACOBIAN_WEIGHT, "{}", quality(&x));
|
||||
}
|
||||
|
||||
/// The analytic gradient matches central differences of the quality away from ties, on
|
||||
/// either branch.
|
||||
#[test]
|
||||
fn quality_gradient_matches_central_differences() {
|
||||
let sliverish = positive([[0.1, 0.0, 0.05], [1.0, 0.2, 0.0], [0.3, 0.9, 0.1], [0.4, 0.35, 0.7]]);
|
||||
let needle = positive([[0.0, 0.0, 0.0], [1.0, 0.1, 0.0], [0.45, 0.85, 0.05], [0.55, 0.3, 6.0]]);
|
||||
for x in [sliverish, needle] {
|
||||
for i in 0..4 {
|
||||
let (q, g) = quality_gradient(&x, i);
|
||||
assert_eq!(q, quality(&x));
|
||||
let h = 1e-7;
|
||||
for c in 0..3 {
|
||||
let (mut a, mut b) = (x, x);
|
||||
a[i][c] += h;
|
||||
b[i][c] -= h;
|
||||
let fd = (quality(&a) - quality(&b)) / (2.0 * h);
|
||||
assert!((fd - g[c]).abs() < 1e-6, "corner {i} axis {c}: {fd} vs {}", g[c]);
|
||||
}
|
||||
}
|
||||
}
|
||||
assert_eq!(quality_gradient(&[[0.0; 3]; 4], 0).1, [0.0; 3]);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn sliver_scores_low() {
|
||||
let x = positive([[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [1.0, 1.0, 0.01]]);
|
||||
|
||||
+93
-46
@@ -31,7 +31,7 @@ pub fn tet_quality(cx: &TetComplex, vs: [u32; 4]) -> Result<f64> {
|
||||
pub fn worst(cx: &TetComplex, tets: &[u32]) -> Result<f64> {
|
||||
let mut q = f64::INFINITY;
|
||||
for &t in tets {
|
||||
q = q.min(tet_quality(cx, cx.tet(t)?)?);
|
||||
q = q.min(cx.slot_quality(t)?);
|
||||
}
|
||||
Ok(q)
|
||||
}
|
||||
@@ -47,13 +47,17 @@ pub const MAX_RING: usize = 10;
|
||||
pub const BOUNDARY_FLIP_COS: f64 = 0.939_692_620_785_908_4;
|
||||
|
||||
/// Klincsek's table: the best worst quality of a triangulation of the chain `v_i … v_j`
|
||||
/// (with the closing edge `v_i v_j`), and the apex `k` that achieves it.
|
||||
fn klincsek(m: usize, q: &mut dyn FnMut(usize, usize, usize) -> Result<f64>) -> Result<(Vec<f64>, Vec<usize>)> {
|
||||
/// (with the closing edge `v_i v_j`), and the apex `k` that achieves it. Only triangulations
|
||||
/// better than `floor` matter: an entry that cannot beat it holds `floor`, and a triangle is
|
||||
/// evaluated only when its two sub-chains beat it, which skips most of the table when the
|
||||
/// removal cannot pay off (the decision and the triangulation chosen are the same as without
|
||||
/// the floor).
|
||||
fn klincsek(m: usize, floor: f64, q: &mut dyn FnMut(usize, usize, usize) -> Result<f64>) -> Result<(Vec<f64>, Vec<usize>)> {
|
||||
let mut best = vec![f64::INFINITY; m * m];
|
||||
let mut arg = vec![0usize; m * m];
|
||||
for i in (0..m.saturating_sub(2)).rev() {
|
||||
for j in i + 2..m {
|
||||
let mut b = f64::NEG_INFINITY;
|
||||
let mut b = floor;
|
||||
let mut a = i + 1;
|
||||
for k in i + 1..j {
|
||||
let sub = best.get(i * m + k).copied().unwrap_or(f64::NEG_INFINITY).min(best.get(k * m + j).copied().unwrap_or(f64::NEG_INFINITY));
|
||||
@@ -115,12 +119,31 @@ pub fn remove_edge_at(cx: &mut TetComplex, t: u32, a: u32, b: u32) -> Result<boo
|
||||
if cx.curve(a, b).is_some() || !vs.contains(&a) || !vs.contains(&b) {
|
||||
return Ok(false);
|
||||
}
|
||||
// A removal that failed fails again until a tetrahedron around the edge changes.
|
||||
if cx.edge_known_to_fail(a, b) {
|
||||
return Ok(false);
|
||||
}
|
||||
let removed = remove_edge_ring(cx, t, a, b)?;
|
||||
if !removed {
|
||||
cx.mark_edge_failed(a, b);
|
||||
}
|
||||
Ok(removed)
|
||||
}
|
||||
|
||||
/// [`remove_edge_at`] past its checks.
|
||||
fn remove_edge_ring(cx: &mut TetComplex, t: u32, a: u32, b: u32) -> Result<bool> {
|
||||
let Some(ring) = cx.edge_ring_at(t, a, b)? else { return Ok(false) };
|
||||
let m = ring.verts.len();
|
||||
if m < 3 || m > MAX_RING + usize::from(!ring.closed) {
|
||||
return Ok(false);
|
||||
}
|
||||
let v = ring.verts.clone();
|
||||
let mut v = ring.verts.clone();
|
||||
if ring.closed {
|
||||
// Start a closed ring at its smallest vertex, so the triangulation chosen does not
|
||||
// depend on the tetrahedron the ring was walked from.
|
||||
let first = v.iter().enumerate().min_by_key(|(_, id)| **id).map_or(0, |(i, _)| i);
|
||||
v.rotate_left(first);
|
||||
}
|
||||
let x = |id: u32| cx.xyz(id);
|
||||
let (xa, xb) = (x(a)?, x(b)?);
|
||||
// An open ring runs from boundary face (a, b, v_0) to (a, b, v_{m−1}): both must lie on one
|
||||
@@ -150,7 +173,7 @@ pub fn remove_edge_at(cx: &mut TetComplex, t: u32, a: u32, b: u32) -> Result<boo
|
||||
Ok(measure::quality(&[pi, pj, pk, xa]).min(measure::quality(&[pi, pk, pj, xb])))
|
||||
};
|
||||
// Triangulate the polygon v_0 … v_{m−1}, closed by the edge v_{m−1} v_0.
|
||||
let (best, arg) = klincsek(m, &mut q)?;
|
||||
let (best, arg) = klincsek(m, before, &mut q)?;
|
||||
let total = best.get(m - 1).copied().unwrap_or(f64::NEG_INFINITY);
|
||||
if !(total > before) {
|
||||
return Ok(false);
|
||||
@@ -199,37 +222,59 @@ pub fn remove_face(cx: &mut TetComplex, t: u32, k: usize) -> Result<bool> {
|
||||
let us = cx.tet(u)?;
|
||||
let face = face_toward(vs, k);
|
||||
let b = us.iter().copied().find(|v| !face.contains(v)).unwrap_or(NONE);
|
||||
if b == NONE || cx.curve(a, b).is_some() || cx.edge_tet(a, b)?.is_some() {
|
||||
if b == NONE || cx.curve(a, b).is_some() || cx.face_known_to_fail(a, b, face) || cx.edge_tet(a, b)?.is_some() {
|
||||
return Ok(false);
|
||||
}
|
||||
let removed = remove_faces_between(cx, t, u, a, b, face)?;
|
||||
if !removed {
|
||||
cx.mark_face_failed(a, b, face);
|
||||
}
|
||||
Ok(removed)
|
||||
}
|
||||
|
||||
/// [`remove_face`] past its checks: tetrahedra `t` (apex `a`) and `u` (apex `b`) share `face`.
|
||||
fn remove_faces_between(cx: &mut TetComplex, t: u32, u: u32, a: u32, b: u32, face: [u32; 3]) -> Result<bool> {
|
||||
// Sandwiched faces between a and b: faces f with tetrahedra (a ∪ f) and (b ∪ f). Grow the
|
||||
// set from the given face across shared edges while every new tetrahedron (a, b, e) on the
|
||||
// ring stays better than the worst removed one; take the best prefix found (a greedy form of
|
||||
// Shewchuk's optimal search, which visits the same candidates).
|
||||
let mut faces: Vec<[u32; 3]> = vec![face];
|
||||
// Each face with the tetrahedron (a ∪ face).
|
||||
let mut faces: Vec<([u32; 3], u32)> = vec![(face, t)];
|
||||
let mut old: Vec<u32> = vec![t, u];
|
||||
let mut best: Option<(f64, Vec<u32>, Vec<[u32; 4]>)> = None;
|
||||
let mut removed_q = worst(cx, &old)?;
|
||||
let mut made: IdMap<(u32, u32), ([u32; 4], f64)> = IdMap::default();
|
||||
for _round in 0..8 {
|
||||
let removed_q = worst(cx, &old)?;
|
||||
// The boundary edges of the face set (edges in exactly one face), oriented as in the face
|
||||
// seen from a: each gives the new tetrahedron (a, b, e) positive.
|
||||
let mut count: IdMap<(u32, u32), (u32, u32, usize)> = IdMap::default();
|
||||
for f in &faces {
|
||||
// seen from a, with the face's tetrahedron on a's side: each gives the new tetrahedron
|
||||
// (a, b, e) positive.
|
||||
let mut count: IdMap<(u32, u32), (u32, u32, u32, usize)> = IdMap::default();
|
||||
for &(f, ta) in &faces {
|
||||
for i in 0..3 {
|
||||
let (p, q) = (f[i], f[(i + 1) % 3]);
|
||||
let e = count.entry(sorted2(p, q)).or_insert((p, q, 0));
|
||||
e.2 += 1;
|
||||
let e = count.entry(sorted2(p, q)).or_insert((p, q, ta, 0));
|
||||
e.3 += 1;
|
||||
}
|
||||
}
|
||||
let mut ring: Vec<(u32, u32)> = count.values().filter(|e| e.2 == 1).map(|e| (e.0, e.1)).collect();
|
||||
let mut ring: Vec<(u32, u32, u32)> = count.values().filter(|e| e.3 == 1).map(|e| (e.0, e.1, e.2)).collect();
|
||||
ring.sort_unstable();
|
||||
let mut new = Vec::with_capacity(ring.len());
|
||||
let mut q = f64::INFINITY;
|
||||
for &(p, r) in &ring {
|
||||
// face (p, r, x) has a on its positive side, so (p, r, a, b)… orient explicitly.
|
||||
let cand = [p, r, b, a];
|
||||
let cand = if cx.orient(cand)? == Sign::Positive { cand } else { [r, p, b, a] };
|
||||
q = q.min(tet_quality(cx, cand)?);
|
||||
for &(p, r, _) in &ring {
|
||||
// A ring edge stays on the ring until a face grows across it: its tetrahedron is
|
||||
// worked out once.
|
||||
let (cand, cq) = match made.get(&(p, r)) {
|
||||
Some(&c) => c,
|
||||
None => {
|
||||
// face (p, r, x) has a on its positive side, so (p, r, a, b)… orient explicitly.
|
||||
let cand = [p, r, b, a];
|
||||
let cand = if cx.orient(cand)? == Sign::Positive { cand } else { [r, p, b, a] };
|
||||
let c = (cand, tet_quality(cx, cand)?);
|
||||
made.insert((p, r), c);
|
||||
c
|
||||
}
|
||||
};
|
||||
q = q.min(cq);
|
||||
new.push(cand);
|
||||
}
|
||||
if q > removed_q && best.as_ref().is_none_or(|(bq, _, _)| q > *bq) {
|
||||
@@ -237,11 +282,12 @@ pub fn remove_face(cx: &mut TetComplex, t: u32, k: usize) -> Result<bool> {
|
||||
}
|
||||
// Grow: a sandwiched neighbour face across a ring edge.
|
||||
let mut grown = false;
|
||||
for &(p, r) in &ring {
|
||||
if let Some((f, ta, tb)) = sandwiched_across(cx, a, b, p, r, &faces)? {
|
||||
faces.push(f);
|
||||
old.push(ta);
|
||||
old.push(tb);
|
||||
for &(p, r, ta) in &ring {
|
||||
if let Some((f, na, nb)) = sandwiched_across(cx, a, b, p, r, ta, &faces)? {
|
||||
faces.push((f, na));
|
||||
old.push(na);
|
||||
old.push(nb);
|
||||
removed_q = removed_q.min(worst(cx, &[na, nb])?);
|
||||
grown = true;
|
||||
break;
|
||||
}
|
||||
@@ -257,27 +303,28 @@ pub fn remove_face(cx: &mut TetComplex, t: u32, k: usize) -> Result<bool> {
|
||||
|
||||
/// A face `(p, r, x)` sandwiched between `a` and `b` (tetrahedra `a ∪ f` and `b ∪ f` both exist)
|
||||
/// next to edge `pr`, not already in `faces`.
|
||||
fn sandwiched_across(cx: &mut TetComplex, a: u32, b: u32, p: u32, r: u32, faces: &[[u32; 3]]) -> Result<Option<([u32; 3], u32, u32)>> {
|
||||
let Some(ring) = cx.edge_ring(p, r)? else { return Ok(None) };
|
||||
// Tetrahedra around pr containing a, and containing b.
|
||||
for &ta in &ring.tets {
|
||||
let va = cx.tet(ta)?;
|
||||
if !va.contains(&a) {
|
||||
continue;
|
||||
}
|
||||
let x = va.iter().copied().find(|&v| v != a && v != p && v != r).unwrap_or(NONE);
|
||||
if x == NONE || x == b {
|
||||
continue;
|
||||
}
|
||||
let f = [p, r, x];
|
||||
if faces.iter().any(|g| sorted3(*g) == sorted3(f)) {
|
||||
continue;
|
||||
}
|
||||
let k = va.iter().position(|&v| v == a).unwrap_or(0);
|
||||
let tb = cx.adj(ta)?.get(k).copied().unwrap_or(NONE);
|
||||
if tb != NONE && cx.tet(tb)?.contains(&b) {
|
||||
return Ok(Some((f, ta, tb)));
|
||||
}
|
||||
fn sandwiched_across(cx: &TetComplex, a: u32, b: u32, p: u32, r: u32, ta: u32, faces: &[([u32; 3], u32)]) -> Result<Option<([u32; 3], u32, u32)>> {
|
||||
// ta = (a, p, r, y) holds the face set's face on edge pr; the other tetrahedron of a around
|
||||
// pr in a face set's neighbourhood lies across its face (a, p, r).
|
||||
let va = cx.tet(ta)?;
|
||||
let Some(ky) = va.iter().position(|&v| v != a && v != p && v != r) else { return Ok(None) };
|
||||
let n = cx.adj(ta)?.get(ky).copied().unwrap_or(NONE);
|
||||
if n == NONE {
|
||||
return Ok(None);
|
||||
}
|
||||
let vn = cx.tet(n)?;
|
||||
let x = vn.iter().copied().find(|&v| v != a && v != p && v != r).unwrap_or(NONE);
|
||||
if x == NONE || x == b || !vn.contains(&a) {
|
||||
return Ok(None);
|
||||
}
|
||||
let f = [p, r, x];
|
||||
if faces.iter().any(|(g, _)| sorted3(*g) == sorted3(f)) {
|
||||
return Ok(None);
|
||||
}
|
||||
let Some(k) = vn.iter().position(|&v| v == a) else { return Ok(None) };
|
||||
let tb = cx.adj(n)?.get(k).copied().unwrap_or(NONE);
|
||||
if tb != NONE && cx.tet(tb)?.contains(&b) {
|
||||
return Ok(Some((f, n, tb)));
|
||||
}
|
||||
Ok(None)
|
||||
}
|
||||
|
||||
@@ -207,6 +207,42 @@ impl LocalStar {
|
||||
measure::quality(&y)
|
||||
}
|
||||
|
||||
/// The quality of tetrahedron `k` with the vertex at `p`, and its gradient there: analytic
|
||||
/// for a valid tetrahedron, central differences with step `delta` for a flat or inverted one
|
||||
/// (whose signed measure has no single smooth branch).
|
||||
fn quality_gradient_at(&self, k: usize, p: P, delta: f64) -> (f64, P) {
|
||||
let Some((x, i)) = self.tets.get(k) else { return (f64::NEG_INFINITY, [0.0; 3]) };
|
||||
let mut y = *x;
|
||||
if let Some(s) = y.get_mut(*i) {
|
||||
*s = p;
|
||||
}
|
||||
let (q, g) = measure::quality_gradient(&y, *i);
|
||||
if q > 0.0 {
|
||||
return (q, g);
|
||||
}
|
||||
let mut g = [0.0; 3];
|
||||
for (c, gc) in g.iter_mut().enumerate() {
|
||||
let (mut a, mut b) = (p, p);
|
||||
a[c] += delta;
|
||||
b[c] -= delta;
|
||||
*gc = (self.quality_at(k, a) - self.quality_at(k, b)) / (2.0 * delta);
|
||||
}
|
||||
(q, g)
|
||||
}
|
||||
|
||||
/// The worst quality at `p` over the tetrahedra in `order`, or the first one found that is
|
||||
/// not above `floor` (then the move is refused whatever the others are).
|
||||
fn worst_above(&self, order: &[usize], p: P, floor: f64) -> f64 {
|
||||
let mut w = f64::INFINITY;
|
||||
for &k in order {
|
||||
w = w.min(self.quality_at(k, p));
|
||||
if !(w > floor) {
|
||||
return w;
|
||||
}
|
||||
}
|
||||
w
|
||||
}
|
||||
|
||||
fn worst_at(&self, p: P) -> f64 {
|
||||
(0..self.tets.len()).map(|k| self.quality_at(k, p)).fold(f64::INFINITY, f64::min)
|
||||
}
|
||||
@@ -215,7 +251,15 @@ impl LocalStar {
|
||||
/// Optimization-based smoothing of `v`: raise the worst quality of its star (see the module
|
||||
/// docs). Returns whether it moved.
|
||||
pub fn smooth_vertex(cx: &mut TetComplex, ctx: &Context, v: u32) -> Result<bool> {
|
||||
if !cx.is_used(v) || cx.kind(v) == VertexKind::Corner {
|
||||
smooth_vertex_below(cx, ctx, v, f64::INFINITY)
|
||||
}
|
||||
|
||||
/// [`smooth_vertex`], only when the worst tetrahedron around `v` is below `good` (Klingner and
|
||||
/// Shewchuk smooth only around tetrahedra below their threshold).
|
||||
pub fn smooth_vertex_below(cx: &mut TetComplex, ctx: &Context, v: u32, good: f64) -> Result<bool> {
|
||||
// The optimization is a function of the star alone: one that found no move finds none
|
||||
// again until a tetrahedron around `v` changes.
|
||||
if !cx.is_used(v) || cx.kind(v) == VertexKind::Corner || cx.smoothing_known_to_fail(v) {
|
||||
return Ok(false);
|
||||
}
|
||||
let star = cx.star(v)?;
|
||||
@@ -237,25 +281,19 @@ pub fn smooth_vertex(cx: &mut TetComplex, ctx: &Context, v: u32) -> Result<bool>
|
||||
if !(len > 0.0 && len.is_finite()) {
|
||||
return Ok(false);
|
||||
}
|
||||
let delta = 1e-6 * len;
|
||||
let mut x = x0;
|
||||
let mut worst = local.worst_at(x);
|
||||
if !(worst < good) {
|
||||
return Ok(false);
|
||||
}
|
||||
let mut moved = false;
|
||||
for _ in 0..8 {
|
||||
// Qualities and gradients (central differences) of the tetrahedra near the worst.
|
||||
// Qualities and analytic gradients of the tetrahedra near the worst.
|
||||
let mut qs = Vec::with_capacity(local.tets.len());
|
||||
for k in 0..local.tets.len() {
|
||||
let q = local.quality_at(k, x);
|
||||
let mut g = [0.0; 3];
|
||||
if q <= worst + 0.1 {
|
||||
for (c, gc) in g.iter_mut().enumerate() {
|
||||
let (mut a, mut b) = (x, x);
|
||||
a[c] += delta;
|
||||
b[c] -= delta;
|
||||
*gc = (local.quality_at(k, a) - local.quality_at(k, b)) / (2.0 * delta);
|
||||
}
|
||||
}
|
||||
qs.push((q, restrict(g, f)));
|
||||
let g = if q <= worst + 0.1 { restrict(local.quality_gradient_at(k, x, 1e-6 * len).1, f) } else { [0.0; 3] };
|
||||
qs.push((q, g));
|
||||
}
|
||||
let active: Vec<P> = qs.iter().filter(|s| s.0 <= worst + 1e-4).map(|s| s.1).collect();
|
||||
let d = min_norm(&active);
|
||||
@@ -272,10 +310,15 @@ pub fn smooth_vertex(cx: &mut TetComplex, ctx: &Context, v: u32) -> Result<bool>
|
||||
alpha = alpha.min((q - worst) / (dd - gd));
|
||||
}
|
||||
}
|
||||
// Trial positions check the worst tetrahedra first and stop at the first that does not
|
||||
// beat the current worst.
|
||||
let mut order: Vec<usize> = (0..qs.len()).collect();
|
||||
let q_of = |i: usize| qs.get(i).map_or(f64::INFINITY, |s| s.0);
|
||||
order.sort_by(|&i, &j| q_of(i).total_cmp(&q_of(j)).then(i.cmp(&j)));
|
||||
let mut accepted = false;
|
||||
for _ in 0..6 {
|
||||
let Some(p) = ctx.place(cx.kind(v), add(x, scale(d, alpha))) else { break };
|
||||
let q = local.worst_at(p.coords());
|
||||
let q = local.worst_above(&order, p.coords(), worst);
|
||||
if p.coords() != x && q > worst && cx.move_is_valid(v, &p, &star)? {
|
||||
cx.move_point(v, p)?;
|
||||
x = p.coords();
|
||||
@@ -290,5 +333,8 @@ pub fn smooth_vertex(cx: &mut TetComplex, ctx: &Context, v: u32) -> Result<bool>
|
||||
break;
|
||||
}
|
||||
}
|
||||
if !moved {
|
||||
cx.mark_smoothing_failed(v);
|
||||
}
|
||||
Ok(moved)
|
||||
}
|
||||
@@ -433,7 +433,17 @@ mod determinism {
|
||||
let (p, t, f) = washer();
|
||||
let m = mesh(&p, &t, &f);
|
||||
assert!(m.elements.len() > 100, "{:?}", m.report);
|
||||
assert_eq!(fingerprint(&m), 6_285_915_972_108_509_003, "{} tets, {} nodes", m.elements.len(), m.points.len());
|
||||
// No needles either: the scaled Jacobian analysts check is above the skill's limit.
|
||||
let sj = m
|
||||
.elements
|
||||
.iter()
|
||||
.map(|e| {
|
||||
let xyz: Vec<[f64; 3]> = e.iter().map(|&v| m.points[v as usize]).collect();
|
||||
crate::quality::element_quality(m.kind, &xyz).unwrap().get(crate::quality::Metric::ScaledJacobian).unwrap()
|
||||
})
|
||||
.fold(f64::INFINITY, f64::min);
|
||||
assert!(sj >= 0.3, "min scaled Jacobian {sj}");
|
||||
assert_eq!(fingerprint(&m), 13_918_536_535_603_991_107, "{} tets, {} nodes", m.elements.len(), m.points.len());
|
||||
}
|
||||
|
||||
/// Listing the triangles, their corners and the vertices in another order, and numbering
|
||||
@@ -497,3 +507,62 @@ fn the_gate_refuses_a_detached_tetrahedron() {
|
||||
assert!(e.contains("share no face") && e.contains("(1, 2)"), "{e}");
|
||||
assert!(super::linear_gate(&p, &tets[..1], v / 2.0).is_ok());
|
||||
}
|
||||
|
||||
/// Undoing an edge removal or a face removal puts every slot back as it was.
|
||||
#[test]
|
||||
fn rollback_restores_every_slot() {
|
||||
use crate::cdt::{CdtOptions, tetrahedralize};
|
||||
use crate::surface::{BoundarySurface, SurfaceOptions};
|
||||
let (p, t, f) = determinism::washer();
|
||||
let s = BoundarySurface::new(&p, &t, &f, &SurfaceOptions::default()).unwrap();
|
||||
let cdt = tetrahedralize(&s, &CdtOptions::default()).unwrap();
|
||||
let (mut cx, _) = crate::tet::complex_from_cdt(&s, &cdt, &crate::tet::QualityOptions::default()).unwrap();
|
||||
let state = |cx: &TetComplex| -> Vec<([u32; 4], [u32; 4], [u32; 4])> {
|
||||
(0..cx.slots() as u32).map(|t| (cx.tet(t).unwrap(), cx.adj(t).unwrap(), cx.tags(t).unwrap())).collect()
|
||||
};
|
||||
let mut changed = 0;
|
||||
for (t, vs) in cx.live_tets().collect::<Vec<_>>() {
|
||||
for &(i, j) in &EDGES {
|
||||
for op in 0..2 {
|
||||
let before = state(&cx);
|
||||
cx.journal.recording = true;
|
||||
let len = cx.journal.steps.len();
|
||||
let did = if op == 0 { remove_edge_at(&mut cx, t, vs[i], vs[j]).unwrap() } else { remove_face(&mut cx, t, i).unwrap() };
|
||||
changed += usize::from(did);
|
||||
cx.rollback(len).unwrap();
|
||||
cx.journal.steps.truncate(len);
|
||||
cx.journal.recording = false;
|
||||
assert_eq!(state(&cx), before, "tet {t} edge {i}{j} op {op}");
|
||||
}
|
||||
}
|
||||
}
|
||||
assert!(changed > 0);
|
||||
}
|
||||
|
||||
/// The cached qualities match a fresh computation after improvement (smoothing, flips,
|
||||
/// insertions and the rollbacks of rejected insertions).
|
||||
#[test]
|
||||
fn cached_qualities_match_fresh_ones() {
|
||||
use crate::cdt::{CdtOptions, tetrahedralize};
|
||||
use crate::size::UniformSize;
|
||||
use crate::surface::{BoundarySurface, SurfaceOptions};
|
||||
use crate::tet::geometry::Tessellation;
|
||||
use crate::tet::refine::Context;
|
||||
let (p, t, f) = determinism::washer();
|
||||
let s = BoundarySurface::new(&p, &t, &f, &SurfaceOptions::default()).unwrap();
|
||||
let cdt = tetrahedralize(&s, &CdtOptions::default()).unwrap();
|
||||
let (mut cx, curve_faces) = crate::tet::complex_from_cdt(&s, &cdt, &crate::tet::QualityOptions::default()).unwrap();
|
||||
let (geometry, size) = (Tessellation::new(&s, Vec::new()), UniformSize(3.0));
|
||||
let ctx = Context { size: &size, geometry: &geometry, curve_faces: &curve_faces };
|
||||
// Fill the cache, then change the mesh.
|
||||
for (t, _) in cx.live_tets().collect::<Vec<_>>() {
|
||||
cx.slot_quality(t).unwrap();
|
||||
}
|
||||
let o = crate::tet::improve::ImproveOptions { max_passes: 3, ..Default::default() };
|
||||
let stats = crate::tet::improve::improve(&mut cx, &ctx, &o).unwrap();
|
||||
assert!(stats.attempted > 0 && stats.smoothed > 0, "{stats:?}");
|
||||
for (t, vs) in cx.live_tets().collect::<Vec<_>>() {
|
||||
assert_eq!(cx.slot_quality(t).unwrap().to_bits(), measure::quality(&cx.coords(vs).unwrap()).to_bits());
|
||||
}
|
||||
cx.validate().unwrap();
|
||||
}
|
||||
@@ -534,11 +534,34 @@ impl Reduction {
|
||||
});
|
||||
// A tie whose secondary node is held already: what is left of the row must be a
|
||||
// constraint of its own, not the difference of two near copies of one.
|
||||
if row.source == Source::Tie
|
||||
&& intended.is_none()
|
||||
&& let Some(&(_, c)) = largest
|
||||
&& c.abs() < NEARLY_REPEATED * gross
|
||||
{
|
||||
let nearly_repeated = row.source == Source::Tie && intended.is_none() && largest.is_some_and(|t| t.1.abs() < NEARLY_REPEATED * gross);
|
||||
// Held by a prescribed value (the secondary surface meeting a boundary condition or an
|
||||
// axis-aligned symmetry plane: a crosspoint), the row can only constrain main degrees
|
||||
// of freedom, and nearly repeated, what it adds is the mortar condition of a node
|
||||
// whose motion is already given: its weights on main nodes off the boundary do not
|
||||
// vanish (the dual basis spreads over the node's faces), so after the substitutions
|
||||
// a small but real constraint is left, not round-off (1.3e-6 of the coefficients that
|
||||
// cancelled to give it where a tet face meets a symmetry plane, #540; round-off would
|
||||
// be ~1e-14 there). Eliminating it divides by that smallness; dropping it leaves the
|
||||
// tie off by that fraction at the crosspoint, below the discretization error. So it
|
||||
// is dropped, with the same threshold that calls it a near repeat. A secondary degree
|
||||
// of freedom another constraint already eliminated (a node secondary in two ties, a
|
||||
// tie and a coupling) is a real conflict and still refused (#389).
|
||||
if nearly_repeated && state.get(row.pivot) == Some(&State::Fixed) {
|
||||
// Still consistent with the prescribed values, or a contradiction.
|
||||
let value = |t: &(usize, f64)| t.1 * values.get(t.0).copied().unwrap_or(0.0);
|
||||
let residual: f64 = terms.iter().filter(|t| state.get(t.0) == Some(&State::Fixed)).map(value).sum();
|
||||
let size: f64 = terms.iter().map(|t| value(t).abs()).sum();
|
||||
if residual.abs() > 1e-9 * size.max(f64::MIN_POSITIVE) && residual.abs() > 0.0 {
|
||||
return Err(Error::Invalid(format!(
|
||||
"{} `{}` contradicts the prescribed values (boundary conditions or couplings hold its nodes in a way it cannot follow)",
|
||||
row.source.what(),
|
||||
row.coupling
|
||||
)));
|
||||
}
|
||||
continue;
|
||||
}
|
||||
if nearly_repeated {
|
||||
let mut others: Vec<String> = row
|
||||
.terms
|
||||
.iter()
|
||||
@@ -858,6 +881,48 @@ impl Reduction {
|
||||
mod tests {
|
||||
use super::*;
|
||||
|
||||
fn tie(terms: Vec<(usize, f64)>, pivot: usize) -> Row<'static> {
|
||||
Row { terms, pivot, coupling: "glue", source: Source::Tie }
|
||||
}
|
||||
|
||||
/// Dofs `held` prescribed, dof 1 at `one` and the others at 0.
|
||||
fn reduce(rows: &[Row], held: &[usize], one: f64) -> Result<Reduction> {
|
||||
let layout = DofLayout::uniform(0..2, DofSet::TRANSLATIONS);
|
||||
let budget = Budget { max_bytes: usize::MAX, pattern_entries: 0, matrices: 1 };
|
||||
Reduction::new(
|
||||
&layout,
|
||||
std::iter::empty(),
|
||||
&[0, 1],
|
||||
rows,
|
||||
|active| Ok((0..active.len()).map(|d| held.contains(&d).then_some(if d == 1 { one } else { 0.0 })).collect()),
|
||||
budget,
|
||||
)
|
||||
}
|
||||
|
||||
/// A tie row of a secondary degree of freedom held by a boundary condition (dof 1), which
|
||||
/// after substituting an earlier tie row (dof 0 following dofs 2 and 3) nearly repeats it, is
|
||||
/// dropped (a crosspoint); dof 3 stays free.
|
||||
#[test]
|
||||
fn a_tie_row_on_a_held_degree_of_freedom_is_dropped() {
|
||||
let rows = [tie(vec![(0, 1.0), (2, -0.5), (3, -0.5)], 0), tie(vec![(1, 1.0), (0, -1.0), (2, 0.5), (3, 0.5 - 1e-7)], 1)];
|
||||
let r = reduce(&rows, &[1, 2], 0.0).unwrap();
|
||||
assert_eq!(r.neq(), 3, "dofs 3, 4, 5 free");
|
||||
}
|
||||
|
||||
/// Nearly repeated but inconsistent: on a secondary degree of freedom another tie already
|
||||
/// eliminated (a node secondary in two ties) it is refused (#389), and on a held one whose
|
||||
/// prescribed values it contradicts too.
|
||||
#[test]
|
||||
fn a_near_repeat_on_an_eliminated_degree_of_freedom_is_refused() {
|
||||
let rows = [tie(vec![(0, 1.0), (2, -0.5), (3, -0.5)], 0), tie(vec![(0, 1.0), (2, -0.5), (3, -0.5 - 1e-7)], 0)];
|
||||
let e = reduce(&rows, &[2], 0.0).err().map(|e| e.to_string()).unwrap_or_default();
|
||||
assert!(e.contains("almost alike"), "{e}");
|
||||
// Dof 1 held at 1 and dof 2 at 0 while the row says dof 1 follows them.
|
||||
let rows = [tie(vec![(0, 1.0), (2, -0.5), (3, -0.5)], 0), tie(vec![(1, 1.0), (0, -1.0), (2, 0.5), (3, 0.5 - 1e-7)], 1)];
|
||||
let e = reduce(&rows, &[1, 2], 1.0).err().map(|e| e.to_string()).unwrap_or_default();
|
||||
assert!(e.contains("contradicts"), "{e}");
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn reduced_matrix_budget_refuses_byte_counts_that_overflow() {
|
||||
for matrices in [1, 2] {
|
||||
|
||||
@@ -52,10 +52,6 @@ const PLAIN_CUP: &str = r#"[
|
||||
/// Why a part fails the skill's limits today, read from its numbers and its section picture.
|
||||
/// Shown only while the part fails.
|
||||
const CAUSES: &[(&str, &str)] = &[
|
||||
(
|
||||
"filleted-bracket",
|
||||
"no structured method takes the whole body (the plan has no cuts for the fillet beside the holes), and tet10 at the plan's size leaves slivers (min SJ 0.04, 1709 marked). The same part partitioned into mapped cells (`filleted-bracket-blocks`) passes: partition (§3.4).",
|
||||
),
|
||||
(
|
||||
"loft-transition",
|
||||
"the lofted sides are swept one to one, and the elements where the section changes shape are skewed (min SJ 0.19). Relaxing the four marked elements did not raise it: the blocks are wrong, not the nodes (§6.3). A loft has no template yet (`unmapped`).",
|
||||
@@ -64,7 +60,6 @@ const CAUSES: &[(&str, &str)] = &[
|
||||
"sharp-wedge",
|
||||
"the wedge's sharp edge is swept as hex, so every element tapers to it (min SJ 0.01, edge ratio 150). A column of wedge6 at the edge (§4) raises the SJ but not the ratio: any element with a corner on a 2.9° edge has sides at least 1/tan 2.9° ≈ 20 times its width there, twice the limit. The edge has to be blunted (defeatured) or the limit waived along it.",
|
||||
),
|
||||
("sphere-on-plate", "no structured method takes the plate with the sunk sphere, and tet10 leaves slivers (min SJ 0.12, 662 marked)."),
|
||||
(
|
||||
"analytic-sphere",
|
||||
"no structured method takes a sphere (it has no face to sweep from or axis to revolve about), so the plan does not apply and the fallback's sphere template meshes it at the plan's size (2.45 mm, which keeps a sphere's chords to 0.75 % on paper), 1.1 % short of the volume: the template's elements at the cube's eight corners span more of the sphere than the size. A sphere template in the plan (§3) would size its caps; hex20 (`order: 2`) closes it too.",
|
||||
|
||||
@@ -398,13 +398,21 @@ Delaunay methods with Voronoi/ODT improvement, then (2) all-hex and hex-dominant
|
||||
- **Implemented** in `fieldcraft-mesh::tet::{ops, improve}` from the two papers: Klincsek's dynamic programming for edge
|
||||
removal (interior edges; boundary edges of one flat B-rep face end in a 2-2 flip), multi-face removal grown greedily
|
||||
over the sandwiched faces, nonsmooth optimization smoothing (active set, smallest vector of the convex hull of the
|
||||
active gradients, line search; gradients by central differences), and the composite vertex insertion with rollback
|
||||
active gradients, line search; analytic gradients of the active measure, as §3.1 and Freitag–Ollivier-Gooch 1997
|
||||
use), and the composite vertex insertion with rollback
|
||||
(barycentres of boundary faces, of the tetrahedron, midpoints of boundary edges, and a point inside off a flat
|
||||
boundary edge). Deviations: the inserted vertex takes the Delaunay cavity rather than the optimal cut of §3.4; each pass
|
||||
smooths and flips the tetrahedra below sin 40° and inserts only when the worst did not improve; passes and the
|
||||
insertion list are capped. Objective: the minimum dihedral sine. Every operation is a checked replacement of a region
|
||||
insertion list are capped. Objective: the smaller of the minimum dihedral sine and the scaled Jacobian at the worst
|
||||
corner weighted so that a Jacobian of 0.3 counts as a 30° dihedral angle: the sine misses needles (every dihedral
|
||||
angle fine, one corner nearly flat; Shewchuk 2002); the scaled Jacobian analysts check (Knupp 2001, B6.6) catches
|
||||
them, and with the weight it decides only near that limit. Every operation is a checked replacement of a region
|
||||
(positive by exact `orient3d`, faces matched with opposite orientation, removed vertices fully inside), so the mesh is
|
||||
valid after each step.
|
||||
valid after each step. Speed (#84): an operation (topological, per tetrahedron, edge and face; or a smoothing) that
|
||||
failed is not tried again until a tetrahedron it reads changes (a change clock per vertex; exact, because the
|
||||
operations are functions of that region and the quality is computed on the corners in a canonical order), live
|
||||
tetrahedra's qualities are cached, Klincsek's table starts at the quality to beat, and the composite insertion runs
|
||||
its flips and neighbour smoothing only where a tetrahedron is below the threshold.
|
||||
|
||||
**A3.7 Cheng, Dey, Edelsbrunner, Facello, Teng 2000, Sliver exudation.** **Verdict: reference.**
|
||||
- **Citation.** J. ACM 47(5):883–904, 2000. DOI [10.1145/355483.355487](https://doi.org/10.1145/355483.355487).
|
||||
|
||||
@@ -31,7 +31,7 @@ Times are wall clock on the machine that ran it.
|
||||
| bolted-assembly | plan (sweep), mesh.hex.generate (auto) | – | 6272 (6272 hex8) | 0.583 | 0.939 | 3.6 | 2.0e-4 | 0 | 0 (not relaxed) | 8 coincident, 0 cracks | – | 3.7 | pass | [x=50.000](skill/bolted-assembly.png) |
|
||||
| cantilever-i-beam | plan (sweep) | – | 61600 (61600 hex8) | 1.000 | 1.000 | 1.0 | 0.0e0 | 0 | 0 (not relaxed) | none | – | 8.5 | pass | [x=100.000](skill/cantilever-i-beam.png) |
|
||||
| cross-cylinders | none | – | – | – | – | – | – | – | – | – | – | 27.2 | no mesh | – |
|
||||
| filleted-bracket | tet10 | – | 21525 (21525 tet10) | 0.021 | 0.535 | 13.5 | -7.0e-6 | 0 | 1492 (not relaxed) | none | – | 15.9 | fail (SJ, ratio) | [x=30.000](skill/filleted-bracket.png) |
|
||||
| filleted-bracket | tet10 | – | 20020 (20020 tet10) | 0.320 | 0.562 | 4.1 | -6.8e-6 | 0 | 0 (not relaxed) | none | – | 12.7 | pass (fallback) | [x=30.000](skill/filleted-bracket.png) |
|
||||
| filleted-bracket-blocks | plan (sweep) | – | 75410 (75410 hex8) | 0.704 | 0.954 | 3.5 | 3.7e-5 | 0 | 0 (not relaxed) | none (after merging 2348 nodes) | – | 6.9 | pass | [x=30.000](skill/filleted-bracket-blocks.png) |
|
||||
| i-beam | plan (sweep) | – | 61600 (61600 hex8) | 1.000 | 1.000 | 1.0 | 0.0e0 | 0 | 0 (not relaxed) | none | – | 4.5 | pass | [z=100.000](skill/i-beam.png) |
|
||||
| loft-transition | plan (sweep) | – | 2873 (2873 hex8) | 0.189 | 0.959 | 1.7 | -2.1e-3 | 0 | 4 → 4 (4 relaxed, 0 undone; min SJ 0.189 → 0.189, max ratio 1.7 → 1.7) | none | – | 0.9 | fail (SJ) | [z=25.000](skill/loft-transition.png) |
|
||||
@@ -40,7 +40,7 @@ Times are wall clock on the machine that ran it.
|
||||
| rounded-box | plan (project) | 2 plane(s) | 11520 (11520 hex8) | 0.717 | 0.973 | 4.0 | -4.1e-3 | 0 | 0 (not relaxed) | none | – | 8.7 | pass | [x=20.000](skill/rounded-box.png) |
|
||||
| shallow-groove | plan (sweep) | – | 1960 (1960 hex8) | 0.991 | 0.999 | 1.1 | 2.7e-4 | 0 | 0 (not relaxed) | none | – | 0.2 | pass | [x=30.000](skill/shallow-groove.png) |
|
||||
| sharp-wedge | plan (sweep) | – | 1674 (1674 hex8) | 0.014 | 0.301 | 149.6 | 2.7e-15 | 0 | 1524 (not relaxed) | none | – | 0.2 | fail (SJ, ratio) | [x=40.000](skill/sharp-wedge.png) |
|
||||
| sphere-on-plate | tet10 | – | 23799 (23799 tet10) | 0.092 | 0.587 | 5.9 | -8.7e-5 | 0 | 679 (not relaxed) | none | – | 13.7 | fail (SJ) | [y=0.000](skill/sphere-on-plate.png) |
|
||||
| sphere-on-plate | tet10 | – | 24175 (24175 tet10) | 0.341 | 0.604 | 4.2 | -9.6e-5 | 0 | 0 (not relaxed) | none | – | 6.5 | pass (fallback) | [y=0.000](skill/sphere-on-plate.png) |
|
||||
| spline-vase | plan (sweep) | – | 4536 (4536 hex8) | 0.665 | 0.945 | 9.9 | -6.6e-3 | 0 | 0 (not relaxed) | none | – | 2.3 | pass | [z=-35.000](skill/spline-vase.png) |
|
||||
| stepped-block | plan (cells) | – | 1416 (1416 hex8) | 1.000 | 1.000 | 1.3 | -4.8e-15 | 0 | 0 (not relaxed) | none | – | 1.6 | pass | [x=30.000](skill/stepped-block.png) |
|
||||
| thin-ribs | plan (sweep) | – | 301308 (301308 hex8) | 0.999 | 1.000 | 1.3 | -3.4e-12 | 0 | 0 (not relaxed) | none | – | 34.7 | pass | [x=50.000](skill/thin-ribs.png) |
|
||||
@@ -75,7 +75,6 @@ limits the cause read from its numbers and section picture.
|
||||
- Bracket: mesh.plan.apply: mesh plan `Plan-Bracket`: the cells ["Bracket"] (cuts []) do not mesh all-hex: body `Bracket`: not supported yet: no structured hex method fits this body (extrude: cylindrical face 8 is not parallel to the direction; sweep: from face 1: lin…
|
||||
- Bracket: mesh.hex.generate: body `Bracket`: not supported yet: no structured hex method fits this body (extrude: cylindrical face 8 is not parallel to the direction; sweep: from face 1: linking face 6 has 2 boundary loops (periodic or holed linking faces are not swept…
|
||||
- Bracket: mesh.hex.project (sphere): body `Bracket`: meshing failed: 91 of 5104 elements are inverted, twisted or degenerate after projection (certified Jacobian; smallest corner scaled Jacobian -0.760): block 1 cell [3, 3, 0] near (12.519, 11.405, 1.604); block 1 cell [3, 4, …
|
||||
- **fail (SJ, ratio)**: no structured method takes the whole body (the plan has no cuts for the fillet beside the holes), and tet10 at the plan's size leaves slivers (min SJ 0.04, 1709 marked). The same part partitioned into mapped cells (`filleted-bracket-blocks`) passes: partition (§3.4).
|
||||
- `loft-transition` (corpus: A 40 mm square lofted to an R15 circle 50 mm above: a ruled free-form surface between unlike sections.):
|
||||
- **fail (SJ)**: the lofted sides are swept one to one, and the elements where the section changes shape are skewed (min SJ 0.19). Relaxing the four marked elements did not raise it: the blocks are wrong, not the nodes (§6.3). A loft has no template yet (`unmapped`).
|
||||
- `sharp-wedge` (corpus: A wedge whose faces meet at under 3 degrees: a sharp edge with a needle-thin cross-section.):
|
||||
@@ -84,7 +83,6 @@ limits the cause read from its numbers and section picture.
|
||||
- Plate: mesh.plan.apply: mesh plan `Plan-Plate`: the cells ["Plate"] (cuts []) do not mesh all-hex: body `Plate`: not supported yet: no structured hex method fits this body (extrude: face 6 is curved across the direction; sweep: from face 0: 6 faces are not on the …
|
||||
- Plate: mesh.hex.generate: body `Plate`: not supported yet: no structured hex method fits this body (extrude: face 6 is curved across the direction; sweep: from face 0: 6 faces are not on the sweep path; from face 5: linking face 11 has 2 boundary loops (periodic or …
|
||||
- Plate: mesh.hex.project (sphere): body `Plate`: meshing failed: the mesh volume 32917.42909007507 is far from the body's 69353.30815403847: the blocks do not cover the part (every outer block face must stand for part of its surface)
|
||||
- **fail (SJ)**: no structured method takes the plate with the sunk sphere, and tet10 leaves slivers (min SJ 0.12, 662 marked).
|
||||
- `analytic-sphere` (corpus: Sphere r 20: curvature everywhere, no edges.):
|
||||
- Sphere: mesh.plan.apply: mesh plan `Plan-Sphere`: the cells ["Sphere"] (cuts []) do not mesh all-hex: body `Sphere`: not supported yet: no structured hex method fits this body (extrude: no direction to try; sweep: no planar face to sweep from; revolve: no axis to t…
|
||||
- Sphere: mesh.hex.generate: body `Sphere`: not supported yet: no structured hex method fits this body (extrude: no direction to try; sweep: no planar face to sweep from; revolve: no axis to try; mapped: meshing failed: meshing failed: mesh validity gate: 8 inverted or…
|
||||
|
||||
Binary file not shown.
|
Before Width: | Height: | Size: 87 KiB After Width: | Height: | Size: 82 KiB |
Binary file not shown.
|
Before Width: | Height: | Size: 102 KiB After Width: | Height: | Size: 105 KiB |
Reference in new issue
Block a user