rustyverse/rustytorch
2
0
Fork 0
Code Issues Pull Requests Actions 2 Packages Projects Releases Wiki Activity

Merge r8h-flat-capsule (R8/R7 phase 2 round 1; default-off, verified)

Browse Source 
Co-Authored-By: Claude Opus 5.5 (1M context) <[email protected]>
This commit is contained in:
Omar Sobh
2026-09-26 04:52:30 -05:00
co-authored by Claude Opus 5.5
parent 2d84dcdbdb 02ab155022
commit a9ba290f03
9 changed files with 401 additions and 69 deletions
Download Patch File Download Diff File Expand all files Collapse all files
crates/specialized/rtx-cfd/src/kernels/cuda/e3_geom.cu
+85 -36
View File
@@ -25,6 +25,8 @@ struct GeomSdf {
int cyl_cut, flag_cut; /* cut to the span */
int npts; /* polyline points (x, y interleaved) */
int nst, ns; /* R8-c plate: stations and points per station (nst 0: the polyline) */
int flat_tip; /* R8-h: 1 = the flat tip (corners rounded to tip_corner); 0 = the capsule */
double tip_corner;
};
struct GeomGrid {
@@ -46,6 +48,24 @@ __device__ __forceinline__ double span_cut(double d2, double z, const GeomSdf& s
return outside + rs_min(rs_max(q1, q2), 0.0) - r;
}
/* R8-h: the tip's axial distance along the last segment a→b from b (plate.rs tip_axial). */
__device__ __forceinline__ double tip_axial(double x, double y, double ax, double ay, double bx, double by)
{
double ex = bx - ax, ey = by - ay;
double l = sqrt(ex * ex + ey * ey);
return ((x - bx) * ex + (y - by) * ey) / l;
}
/* R8-h: the flat tip's rounded cap (plate.rs flat_cap). */
__device__ __forceinline__ double flat_cap(double dl, double a, double r)
{
double q1 = dl + r;
double q2 = a + r;
double m1 = rs_max(q1, 0.0), m2 = rs_max(q2, 0.0);
double outside = sqrt(m1 * m1 + m2 * m2);
return outside + rs_min(rs_max(q1, q2), 0.0) - r;
}
/*
* R8-c: the deformed plate (`plate.rs`, the host twin expression for
* expression): `poly` holds the stations' rows (x, y interleaved, row-major
@@ -57,7 +77,7 @@ __device__ __forceinline__ double span_cut(double d2, double z, const GeomSdf& s
*/
__device__ double plate_dist(double x, double y, double z, const GeomSdf& s,
const double* __restrict__ P, const double* __restrict__ V,
double* vx, double* vy)
double* vx, double* vy, double* axial)
{
int nst = s.nst, ns = s.ns;
const double* Z = P + 2 * (long long) nst * ns;
@@ -97,14 +117,14 @@ __device__ double plate_dist(double x, double y, double z, const GeomSdf& s,
double l2 = ex * ex + ey * ey;
double u = ((x - ax) * ex + (y - ay) * ey) / l2;
if (u < 0.0) u = 0.0;
if (u > 1.0) u = 1.0;
if (u > 1.0 && !(s.flat_tip && m + 2 == ns)) u = 1.0;
double px = ax + u * ex, py = ay + u * ey;
double qx = x - px, qy = y - py;
double d = sqrt(qx * qx + qy * qy);
if (d < best) {
best = d;
mb = m;
ub = u;
ub = u > 1.0 ? 1.0 : u;
qbx = qx;
qby = qy;
}
@@ -118,6 +138,20 @@ __device__ double plate_dist(double x, double y, double z, const GeomSdf& s,
double qn = (qbx * cx + qby * cy) / best;
best = best / sqrt(1.0 + qn * qn);
}
if (s.flat_tip && axial) {
int m = ns - 2;
double ax, ay, bx, by;
if (interp) {
ax = R0[2 * m] + w * (R1[2 * m] - R0[2 * m]);
ay = R0[2 * m + 1] + w * (R1[2 * m + 1] - R0[2 * m + 1]);
bx = R0[2 * m + 2] + w * (R1[2 * m + 2] - R0[2 * m + 2]);
by = R0[2 * m + 3] + w * (R1[2 * m + 3] - R0[2 * m + 3]);
} else {
ax = R0[2 * m]; ay = R0[2 * m + 1];
bx = R0[2 * m + 2]; by = R0[2 * m + 3];
}
*axial = tip_axial(x, y, ax, ay, bx, by);
}
if (V) {
const double* V0 = V + 2 * (long long) k * ns;
const double* V1 = interp ? V0 + 2 * ns : V0;
@@ -144,22 +178,29 @@ __device__ double geom_phi_at(double x, double y, double z, const GeomSdf& s, co
double dc = sqrt(ex0 * ex0 + ey0 * ey0) - s.rc;
if (s.cyl_cut) dc = span_cut(dc, z, s);
/* the capsule: distance to the polyline (or the plate, R8-c) */
double best = 1.0 / 0.0;
if (s.nst > 0) best = plate_dist(x, y, z, s, poly, nullptr, nullptr, nullptr);
else for (int m = 0; m + 1 < s.npts; ++m) {
double ax = poly[2 * m], ay = poly[2 * m + 1];
double bx = poly[2 * m + 2], by = poly[2 * m + 3];
double ex = bx - ax, ey = by - ay;
double l2 = ex * ex + ey * ey;
double u = ((x - ax) * ex + (y - ay) * ey) / l2;
if (u < 0.0) u = 0.0;
if (u > 1.0) u = 1.0;
double px = ax + u * ex, py = ay + u * ey;
double qx = x - px, qy = y - py;
double d = sqrt(qx * qx + qy * qy);
if (d < best) best = d;
double best = 1.0 / 0.0, axial = 0.0;
if (s.nst > 0) best = plate_dist(x, y, z, s, poly, nullptr, nullptr, nullptr, &axial);
else {
for (int m = 0; m + 1 < s.npts; ++m) {
double ax = poly[2 * m], ay = poly[2 * m + 1];
double bx = poly[2 * m + 2], by = poly[2 * m + 3];
double ex = bx - ax, ey = by - ay;
double l2 = ex * ex + ey * ey;
double u = ((x - ax) * ex + (y - ay) * ey) / l2;
if (u < 0.0) u = 0.0;
if (u > 1.0 && !(s.flat_tip && m + 2 == s.npts)) u = 1.0;
double px = ax + u * ex, py = ay + u * ey;
double qx = x - px, qy = y - py;
double d = sqrt(qx * qx + qy * qy);
if (d < best) best = d;
}
if (s.flat_tip) {
int m = s.npts - 2;
axial = tip_axial(x, y, poly[2 * m], poly[2 * m + 1], poly[2 * m + 2], poly[2 * m + 3]);
}
}
double df = best - s.half;
if (s.flat_tip) df = flat_cap(df, axial, s.tip_corner);
if (s.flag_cut) df = span_cut(df, z, s);
double r = s.fillet;
if (r > 0.0 && dc < r && df < r) {
@@ -391,28 +432,36 @@ extern "C" __global__ void e3_geom_gather(
__device__ double body_velocity(double x, double y, double z, int c, const GeomSdf& s,
const double* __restrict__ poly, const double* __restrict__ vel)
{
double best = 1.0 / 0.0, vx = 0.0, vy = 0.0;
if (s.nst > 0) best = plate_dist(x, y, z, s, poly, vel, &vx, &vy);
else for (int m = 0; m + 1 < s.npts; ++m) {
double ax = poly[2 * m], ay = poly[2 * m + 1];
double bx = poly[2 * m + 2], by = poly[2 * m + 3];
double ex = bx - ax, ey = by - ay;
double l2 = ex * ex + ey * ey;
double u = ((x - ax) * ex + (y - ay) * ey) / l2;
if (u < 0.0) u = 0.0;
if (u > 1.0) u = 1.0;
double px = ax + u * ex, py = ay + u * ey;
double qx = x - px, qy = y - py;
double d = sqrt(qx * qx + qy * qy);
if (d < best) {
best = d;
double avx = vel[2 * m], avy = vel[2 * m + 1];
double bvx = vel[2 * m + 2], bvy = vel[2 * m + 3];
vx = avx + u * (bvx - avx);
vy = avy + u * (bvy - avy);
double best = 1.0 / 0.0, vx = 0.0, vy = 0.0, axial = 0.0;
if (s.nst > 0) best = plate_dist(x, y, z, s, poly, vel, &vx, &vy, &axial);
else {
for (int m = 0; m + 1 < s.npts; ++m) {
double ax = poly[2 * m], ay = poly[2 * m + 1];
double bx = poly[2 * m + 2], by = poly[2 * m + 3];
double ex = bx - ax, ey = by - ay;
double l2 = ex * ex + ey * ey;
double u = ((x - ax) * ex + (y - ay) * ey) / l2;
if (u < 0.0) u = 0.0;
if (u > 1.0 && !(s.flat_tip && m + 2 == s.npts)) u = 1.0;
double px = ax + u * ex, py = ay + u * ey;
double qx = x - px, qy = y - py;
double d = sqrt(qx * qx + qy * qy);
if (d < best) {
best = d;
double uv = u > 1.0 ? 1.0 : u;
double avx = vel[2 * m], avy = vel[2 * m + 1];
double bvx = vel[2 * m + 2], bvy = vel[2 * m + 3];
vx = avx + uv * (bvx - avx);
vy = avy + uv * (bvy - avy);
}
}
if (s.flat_tip) {
int m = s.npts - 2;
axial = tip_axial(x, y, poly[2 * m], poly[2 * m + 1], poly[2 * m + 2], poly[2 * m + 3]);
}
}
double df = best - s.half;
if (s.flat_tip) df = flat_cap(df, axial, s.tip_corner);
if (s.flag_cut) df = span_cut(df, z, s);
double ex0 = x - s.cx, ey0 = y - s.cy;
double dc = sqrt(ex0 * ex0 + ey0 * ey0) - s.rc;
crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/body.rs
+6
View File
@@ -31,6 +31,12 @@ pub struct DeviceSdf {
pub half: f64,
/// The root fillet radius (0 = the plain union).
pub fillet: f64,
/// R8-h: the tip. `None` = the capsule (the semicircular tip around the
/// last point, as before); `Some(r_c)` = a FLAT tip: the flag ends in
/// the plane through the last point normal to the last segment, its two
/// corners rounded to radius `r_c` (0 ≤ r_c ≤ `half`; `r_c = half` is
/// the capsule's semicircle pulled back by `half`, geometrically).
pub tip_corner: Option<f64>,
/// The centreline polyline's points (x, y) at `t`.
pub poly: Vec<[f64; 2]>,
/// R6-2 step 2: the centreline's velocity (vx, vy) per point at `t`, the
crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/interface.rs
+1
View File
@@ -708,6 +708,7 @@ mod tests {
r_edge: 0.0,
half,
fillet: 0.0,
tip_corner: None,
poly: Vec::new(),
vel: Vec::new(),
plate: Some(surf),
crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/plate.rs
+65 -12
View File
@@ -131,12 +131,16 @@ fn span_cut(d2: f64, z: f64, s: &DeviceSdf) -> f64 {
}
/// The closest segment of a polyline given by `pt(m)`, `m < n`: the
/// in-plane distance, the segment, its parameter and the offset `(qx, qy)`.
/// in-plane distance, the segment, its parameter (clamped to the segment)
/// and the offset `(qx, qy)`. R8-h: with `ray_last` the last segment is a
/// RAY beyond its end point (its parameter unclamped above for the
/// distance) — the lateral distance the flat tip's cap needs.
#[inline]
fn closest<P: Fn(usize) -> [f64; 2]>(
x: f64,
y: f64,
n: usize,
ray_last: bool,
pt: P,
) -> (f64, usize, f64, [f64; 2]) {
let mut best = f64::INFINITY;
@@ -151,7 +155,7 @@ fn closest<P: Fn(usize) -> [f64; 2]>(
if u < 0.0 {
u = 0.0;
}
if u > 1.0 {
if u > 1.0 && !(ray_last && m + 2 == n) {
u = 1.0;
}
let px = ax + u * ex;
@@ -162,25 +166,63 @@ fn closest<P: Fn(usize) -> [f64; 2]>(
if d < best {
best = d;
mb = m;
ub = u;
ub = if u > 1.0 { 1.0 } else { u };
qb = [qx, qy];
}
}
(best, mb, ub, qb)
}
/// R8-h: the signed distance along the last segment's unit tangent from
/// its end point `b` (the flat tip's plane; > 0 beyond the tip).
#[inline]
fn tip_axial(x: f64, y: f64, a: [f64; 2], b: [f64; 2]) -> f64 {
let ex = b[0] - a[0];
let ey = b[1] - a[1];
let l = (ex * ex + ey * ey).sqrt();
((x - b[0]) * ex + (y - b[1]) * ey) / l
}
/// R8-h: the flat tip — the strip of lateral distance `dl` (distance to
/// the centreline with its last segment a ray, minus the half-thickness)
/// cut by the tip's plane (axial distance `a`) with the two corners rounded
/// to radius `r` (the span cut's rounded intersection). Exact near the tip
/// where the centreline is straight over the corner (the last segment).
#[inline]
fn flat_cap(dl: f64, a: f64, r: f64) -> f64 {
let q1 = dl + r;
let q2 = a + r;
let m1 = rs_max(q1, 0.0);
let m2 = rs_max(q2, 0.0);
let outside = (m1 * m1 + m2 * m2).sqrt();
outside + rs_min(rs_max(q1, q2), 0.0) - r
}
/// The plate's in-plane closest point at `(x, y, z)`, the slope-corrected
/// distance to the mid-surface, and the interpolated velocity there
/// (zero without velocities).
fn plate_closest(p: &PlateSurface, x: f64, y: f64, z: f64) -> (f64, [f64; 2]) {
///
/// R8-h: with `flat` the last segment is a ray (see [`closest`]) and the
/// third value is the tip's axial distance on the interpolated polyline
/// (0 otherwise).
fn plate_closest(p: &PlateSurface, x: f64, y: f64, z: f64, flat: bool) -> (f64, [f64; 2], f64) {
let ns = p.ns;
let (k, w) = p.bracket(z);
let row = |k: usize, m: usize| p.xy[k * ns + m];
let lerp2 =
|a: [f64; 2], b: [f64; 2], w: f64| [a[0] + w * (b[0] - a[0]), a[1] + w * (b[1] - a[1])];
let (best, m, u, q) = match w {
None => closest(x, y, ns, |m| row(k, m)),
Some(w) => closest(x, y, ns, |m| lerp2(row(k, m), row(k + 1, m), w)),
None => closest(x, y, ns, flat, |m| row(k, m)),
Some(w) => closest(x, y, ns, flat, |m| lerp2(row(k, m), row(k + 1, m), w)),
};
let axial = if flat {
let at = |m: usize| match w {
None => row(k, m),
Some(w) => lerp2(row(k, m), row(k + 1, m), w),
};
tip_axial(x, y, at(ns - 2), at(ns - 1))
} else {
0.0
};
let mut d = best;
if let Some(_w) = w {
@@ -211,7 +253,7 @@ fn plate_closest(p: &PlateSurface, x: f64, y: f64, z: f64) -> (f64, [f64; 2]) {
let (a, b) = (at(m), at(m + 1));
[a[0] + u * (b[0] - a[0]), a[1] + u * (b[1] - a[1])]
};
(d, v)
(d, v, axial)
}
impl DeviceSdf {
@@ -231,20 +273,30 @@ impl DeviceSdf {
/// and its surface velocity at the closest point.
#[must_use]
pub fn flag_distance_host(&self, x: f64, y: f64, z: f64) -> (f64, [f64; 2]) {
let (best, v) = match self.plate.as_ref() {
Some(p) => plate_closest(p, x, y, z),
let flat = self.tip_corner.is_some();
let (best, v, axial) = match self.plate.as_ref() {
Some(p) => plate_closest(p, x, y, z, flat),
None => {
let (best, m, u, _) = closest(x, y, self.poly.len(), |m| self.poly[m]);
let n = self.poly.len();
let (best, m, u, _) = closest(x, y, n, flat, |m| self.poly[m]);
let v = if self.vel.len() == self.poly.len() && !self.vel.is_empty() {
let (a, b) = (self.vel[m], self.vel[m + 1]);
[a[0] + u * (b[0] - a[0]), a[1] + u * (b[1] - a[1])]
} else {
[0.0, 0.0]
};
(best, v)
let axial = if flat {
tip_axial(x, y, self.poly[n - 2], self.poly[n - 1])
} else {
0.0
};
(best, v, axial)
}
};
let mut df = best - self.half;
if let Some(rc) = self.tip_corner {
df = flat_cap(df, axial, rc);
}
if self.flag_cut {
df = span_cut(df, z, self);
}
@@ -344,6 +396,7 @@ mod tests {
r_edge: 0.0066,
half: 0.01,
fillet: 0.0,
tip_corner: None,
poly,
vel,
plate,
@@ -406,7 +459,7 @@ mod tests {
vel: Vec::new(),
};
for &(y, zz) in &[(0.53, 0.05), (0.47, 0.15), (0.6, 0.12)] {
let (d, _) = plate_closest(&p, 0.1, y, zz);
let (d, _, _) = plate_closest(&p, 0.1, y, zz, false);
let exact = (y - (0.5 + a * (zz - 0.1))).abs() / (1.0 + a * a).sqrt();
assert!((d - exact).abs() < 1e-15, "{d} vs {exact}");
}
crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/step/device/geom.rs
+5
View File
@@ -75,6 +75,9 @@ struct GeomSdf {
npts: i32,
nst: i32,
ns: i32,
/// R8-h: the flat tip (1) and its corner radius.
flat_tip: i32,
tip_corner: f64,
}
unsafe impl DeviceRepr for GeomSdf {}
unsafe impl ValidAsZeroBits for GeomSdf {}
@@ -151,6 +154,8 @@ fn geom_sdf(sdf: &DeviceSdf) -> GeomSdf {
npts: sdf.poly.len() as i32,
nst: sdf.plate.as_ref().map_or(0, |p| p.z.len() as i32),
ns: sdf.plate.as_ref().map_or(0, |p| p.ns as i32),
flat_tip: i32::from(sdf.tip_corner.is_some()),
tip_corner: sdf.tip_corner.unwrap_or(0.0),
}
}
crates/specialized/rtx-cfd/tests/embedded3_flag_wake.rs
+34 -3
View File
@@ -96,8 +96,35 @@ fn root_fillet() -> f64 {
/// record before this date had the apex 10 mm beyond A (`=0` restores
/// them): at ny 62 on the recorded motion that was drag 253.3 → 240.0 and
/// the lift swing 1,005 → 836 (the overset's 867).
///
/// R8-h: with the flat tip the default inset is 0 (the flat face through A).
fn tip_inset() -> f64 {
env_f("RTX_E3_FLAG_TIP_INSET", FLAG_HALF)
env_f(
"RTX_E3_FLAG_TIP_INSET",
if flat_tip().is_some() { 0.0 } else { FLAG_HALF },
)
}
/// R8-h (2026-09-25): the tip's shape. `RTX_E3_FLAG_TIP=flat` gives the
/// flag a FLAT tip through the centreline's last point (A, the inset
/// defaulting to 0), its corners rounded to `RTX_E3_FLAG_TIP_CORNER`
/// metres (default 0.00125, the 2D overset's recommended line; at most
/// `FLAG_HALF`); unset or `capsule` = the capsule (the semicircular tip).
/// The flat tip's host φ and surface velocity are the device form's
/// (`DeviceSdf::phi_host`, the kernel's arithmetic).
fn flat_tip() -> Option<f64> {
match std::env::var("RTX_E3_FLAG_TIP").as_deref() {
Err(_) | Ok("capsule") => None,
Ok("flat") => {
let rc = env_f("RTX_E3_FLAG_TIP_CORNER", 0.00125);
assert!(
(0.0..=FLAG_HALF).contains(&rc),
"RTX_E3_FLAG_TIP_CORNER {rc} outside [0, {FLAG_HALF}]"
);
Some(rc)
}
Ok(v) => panic!("RTX_E3_FLAG_TIP={v}: flat or capsule"),
}
}
/// Smooth union with a concave fillet of radius `r` (the plain `min` at r = 0).
@@ -491,6 +518,7 @@ fn flag_wake_on_the_device() {
r_edge,
half: FLAG_HALF,
fillet: r_fillet,
tip_corner: flat_tip(),
poly: match (&plate, recorded()) {
(Some(_), _) => Vec::new(),
(None, Some(rec)) => recorded_polyline(rec, t)
@@ -527,9 +555,9 @@ fn flag_wake_on_the_device() {
c.1.clone().expect("sdf")
})
};
let body = if plate_body() {
let body = if plate_body() || flat_tip().is_some() {
// R8-c: the host φ and surface velocity ARE the device form's (the
// kernel's arithmetic on the host).
// kernel's arithmetic on the host); R8-h: the flat tip too.
Body::from_sdf(move |x, y, z, t| sdf_at(t).phi_host(x, y, z))
.with_surface_velocity(move |x, y, z, t| sdf_at(t).velocity_host(x, y, z))
} else {
@@ -587,6 +615,9 @@ fn flag_wake_on_the_device() {
g.cells(),
(t_end / dt).ceil() as usize
);
if let Some(rc) = flat_tip() {
println!(" R8-h: FLAT tip through the centreline's last point, corner radius {rc:.5} m");
}
unsafe { std::env::set_var("RTX_PROFILE", "1") };
let mut device = DeviceStep::new(solver, g);
device.upload(&field);
crates/specialized/rtx-cfd/tests/embedded3_flat_tip.rs
+138
View File
@@ -0,0 +1,138 @@
//! R8-h (2026-09-25): the flat-tipped flag's geometry (`DeviceSdf::tip_corner`,
//! host twin `DeviceSdf::phi_host` = the kernel's `geom_phi_at`).
//!
//! G2: a straight flag (a five-point polyline, root a semicircle, tip flat
//! with corners rounded to `r_c`) at two orientations, cut on the cubic grid
//! at three resolutions: the solid volume and the wall area of the cut
//! geometry against the analytic rounded rectangle; the tip plane through
//! the last point; `r_c = half` = the capsule pulled back by `half`.
use rtx_cfd::solvers::incompressible::embedded3::{Body, CutGeometry, DeviceSdf, Grid};
const HALF: f64 = 0.01;
const LEN: f64 = 0.35;
const P0: [f64; 2] = [0.05, 0.08];
fn flag(theta: f64, tip: Option<f64>, pull_back: f64) -> DeviceSdf {
let (c, s) = (theta.cos(), theta.sin());
let len = LEN - pull_back;
let poly = (0..5)
.map(|m| {
let f = len * m as f64 / 4.0;
[P0[0] + f * c, P0[1] + f * s]
})
.collect();
DeviceSdf {
// the circle far outside the domain
cyl: [-10.0, -10.0, 0.01],
cyl_cut: false,
flag_cut: false,
zc: 0.0,
span: 1.0,
r_edge: 0.0,
half: HALF,
fillet: 0.0,
tip_corner: tip,
poly,
vel: Vec::new(),
plate: None,
}
}
/// The analytic cross-section: the strip `LEN × 2 HALF`, the root's half
/// disk, less the two rounded corners' `(1 − π/4) r²` (capsule: a second
/// half disk instead of the flat end).
fn analytic(tip: Option<f64>) -> (f64, f64) {
let pi = std::f64::consts::PI;
match tip {
Some(r) => (
LEN * 2.0 * HALF + 0.5 * pi * HALF * HALF - 2.0 * (1.0 - 0.25 * pi) * r * r,
2.0 * (LEN - r) + (2.0 * HALF - 2.0 * r) + pi * r + pi * HALF,
),
None => (
LEN * 2.0 * HALF + pi * HALF * HALF,
2.0 * LEN + 2.0 * pi * HALF,
),
}
}
fn cut(sdf: DeviceSdf, h: f64) -> (f64, f64) {
let (nx, ny, nz) = ((0.46 / h).ceil() as usize, (0.26 / h).ceil() as usize, 2);
let g = Grid::cubic(nx, ny, nz, h);
let body = Body::from_sdf(move |x, y, z, _t| sdf.phi_host(x, y, z));
let geo = CutGeometry::build(&body, g, 0.0);
let depth = nz as f64 * h;
let solid = (nx * ny * nz) as f64 * h * h * h - geo.fluid_volume();
let (wall, _) = geo.wall_area_and_closure();
(solid / depth, wall / depth)
}
#[test]
fn flat_tip_plane_passes_through_the_last_point() {
for theta in [0.0, 0.3] {
for rc in [0.0, 0.00125, 0.005, HALF] {
let s = flag(theta, Some(rc), 0.0);
let tip = *s.poly.last().unwrap();
let (c, sn) = (f64::cos(theta), f64::sin(theta));
// the face centre on the surface; 1 mm beyond it at 1 mm, inside by 1 mm
assert!(s.phi_host(tip[0], tip[1], 0.0).abs() < 1e-15);
if rc < HALF {
let beyond = s.phi_host(tip[0] + 1e-3 * c, tip[1] + 1e-3 * sn, 0.0);
assert!((beyond - 1e-3).abs() < 1e-12, "{beyond}");
}
let inside = s.phi_host(tip[0] - 1e-3 * c, tip[1] - 1e-3 * sn, 0.0);
assert!(inside < 0.0);
// the corner's apex region: the rounded corner's arc
let (px, py) = (
tip[0] - rc * c - (HALF - rc) * sn,
tip[1] - rc * sn + (HALF - rc) * c,
);
let d = 2e-3;
let (dx, dy) = ((c - sn) / 2f64.sqrt(), (sn + c) / 2f64.sqrt());
let got = s.phi_host(px + (rc + d) * dx, py + (rc + d) * dy, 0.0);
assert!((got - d).abs() < 1e-12, "corner {rc}: {got} vs {d}");
}
}
}
#[test]
fn flat_tip_at_half_is_the_capsule_pulled_back() {
let a = flag(0.2, Some(HALF), 0.0);
let b = flag(0.2, None, HALF);
let mut worst = 0.0f64;
for i in 0..400 {
for j in 0..200 {
let (x, y) = (0.30 + 0.0003 * i as f64, 0.10 + 0.0010 * j as f64);
let (pa, pb) = (a.phi_host(x, y, 0.0), b.phi_host(x, y, 0.0));
if pa.abs().min(pb.abs()) < 0.02 {
worst = worst.max((pa - pb).abs());
}
}
}
println!(" r_c = half vs capsule pulled back by half: max |Δφ| {worst:.2e}");
assert!(worst < 1e-12, "{worst}");
}
#[test]
fn flat_tip_cut_volume_and_area_converge_to_the_rounded_rectangle() {
println!(" theta tip h solid err wall err");
for theta in [0.0, 0.3] {
for tip in [None, Some(0.0), Some(0.00125), Some(0.005)] {
let (va, pa) = analytic(tip);
let mut errs = Vec::new();
for ny in [62.0, 124.0, 248.0] {
let h = 0.41 / ny;
let (v, p) = cut(flag(theta, tip, 0.0), h);
let (ev, ep) = ((v - va) / va, (p - pa) / pa);
println!(
" {theta:.1} {:>9} {h:.3e} {ev:+.3e} {ep:+.3e}",
tip.map_or("capsule".to_string(), |r| format!("{r:.5}"))
);
assert!(ev.abs() < 1e-2 && ep.abs() < 2e-2, "{ev} {ep} at h {h}");
errs.push((ev, ep));
}
// the volume converges (the rung ny 248 below ny 62's error)
assert!(errs[2].0.abs() < errs[0].0.abs(), "volume {errs:?}");
}
}
}