#![allow(
clippy::disallowed_macros,
reason = "raw-helper: this module calls the platform f32/f64 max/min via crate::raw_minmax::raw_minmax_here!() at sites proven safe by nonzero-const:/both-abs:/int-cast: reasoning, see scripts/raw_minmax_manifest.toml"
)]
crate::raw_minmax::raw_minmax_here!(max_f32, min_f32);
use alice_sdf::interval::{eval_lipschitz, Interval, Vec3Interval};
use alice_sdf::SdfNode;
use glam::Vec3;
#[derive(Debug, Clone, Copy, PartialEq)]
pub enum Priority {
Hard,
Soft(f32),
}
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
#[non_exhaustive]
pub enum Evidence {
Proved,
Witnessed,
Modelled {
model: &'static str,
},
}
impl Evidence {
#[must_use]
pub const fn is_proof(self) -> bool {
!matches!(self, Self::Modelled { .. })
}
#[must_use]
pub const fn model(self) -> Option<&'static str> {
match self {
Self::Modelled { model } => Some(model),
_ => None,
}
}
}
#[derive(Debug, Clone, PartialEq, Eq)]
pub struct NotProvable {
pub constraint: &'static str,
pub model: &'static str,
}
impl std::fmt::Display for NotProvable {
fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
write!(
f,
"{} は証明も反例も返せない (根拠: {}) ので Priority::Hard を名乗れない — Law::soft で重みを付けて使う",
self.constraint, self.model
)
}
}
impl std::error::Error for NotProvable {}
#[derive(Debug, Clone)]
pub enum Constraint {
NonOverlap {
a: SdfNode,
b: SdfNode,
},
Containment {
inner: SdfNode,
outer: SdfNode,
},
MinThickness {
node: SdfNode,
min_thickness: f32,
},
Stress {
node: SdfNode,
load_points: Vec<(Vec3, f32)>,
min_thickness_factor: f32,
},
Thermal {
node: SdfNode,
heat_sources: Vec<Vec3>,
search_radius: f32,
min_surface_ratio: f32,
},
Contact {
a: SdfNode,
b: SdfNode,
min_distance: f32,
max_distance: f32,
},
Continuity {
node: SdfNode,
seed_point: Vec3,
},
GradientBound {
node: SdfNode,
max_gradient: f32,
probe: f32,
},
Reachable {
node: SdfNode,
from: Vec3,
to: Vec3,
},
VolumeConservation {
before: SdfNode,
after: SdfNode,
relative_tolerance: f32,
},
#[cfg(feature = "physics")]
ThermalField {
node: SdfNode,
sources: Vec<(Vec3, f32)>,
material: alice_physics::transient_thermal::ThermalMaterial,
metres_per_unit: f32,
ambient_c: f32,
max_temperature_c: f32,
duration_s: f32,
},
}
impl Constraint {
#[must_use]
pub const fn name(&self) -> &'static str {
match self {
Self::NonOverlap { .. } => "NonOverlap",
Self::Containment { .. } => "Containment",
Self::MinThickness { .. } => "MinThickness",
Self::Stress { .. } => "Stress",
Self::Thermal { .. } => "Thermal",
Self::Contact { .. } => "Contact",
Self::Continuity { .. } => "Continuity",
Self::GradientBound { .. } => "GradientBound",
Self::Reachable { .. } => "Reachable",
Self::VolumeConservation { .. } => "VolumeConservation",
#[cfg(feature = "physics")]
Self::ThermalField { .. } => "ThermalField",
}
}
#[must_use]
pub const fn evidence_class(&self) -> Evidence {
match self {
Self::Reachable { .. } => Evidence::Proved,
Self::NonOverlap { .. }
| Self::Containment { .. }
| Self::MinThickness { .. }
| Self::Contact { .. }
| Self::GradientBound { .. } => Evidence::Witnessed,
Self::Stress { .. } => Evidence::Modelled {
model: MODEL_STRESS,
},
Self::Thermal { .. } => Evidence::Modelled {
model: MODEL_THERMAL,
},
Self::Continuity { .. } => Evidence::Modelled {
model: MODEL_CONTINUITY,
},
Self::VolumeConservation { .. } => Evidence::Modelled {
model: MODEL_VOLUME,
},
#[cfg(feature = "physics")]
Self::ThermalField { .. } => Evidence::Modelled {
model: MODEL_THERMAL_FIELD,
},
}
}
}
const MODEL_STRESS: &str =
"肉厚 ≥ force × min_thickness_factor の幾何 proxy (材料 / 断面係数 / 降伏応力を持たない)";
const MODEL_THERMAL: &str =
"熱源近傍の 表面近傍セル数 / 内部セル数 比 (|f| を距離として流用、格子解像度依存)";
const MODEL_CONTINUITY: &str =
"格子セル中心の点標本 + 6-connected flood fill (格子より細い接続 / 分離を取りこぼす)";
#[cfg(feature = "physics")]
const MODEL_THERMAL_FIELD: &str = "alice_physics の陽解法 3D 熱伝導 (格子離散化、\
境界条件 h と形状を両端で挟む: 上界 = 断熱 × 内部確定セルのみ、下界 = 等温 × 同)";
const MODEL_VOLUME: &str =
"格子セル中心の点標本による内部セル数 count (格子より細い体積差を取りこぼす)";
#[derive(Debug, Clone)]
pub struct Law {
pub name: String,
pub priority: Priority,
pub constraint: Constraint,
}
impl Law {
pub fn hard(name: impl Into<String>, constraint: Constraint) -> Result<Self, NotProvable> {
if let Some(model) = constraint.evidence_class().model() {
return Err(NotProvable {
constraint: constraint.name(),
model,
});
}
Ok(Self::hard_unchecked(name, constraint))
}
fn hard_unchecked(name: impl Into<String>, constraint: Constraint) -> Self {
debug_assert!(
constraint.evidence_class().is_proof(),
"hard_unchecked にモデル推定の制約 {} が渡された",
constraint.name()
);
Self {
name: name.into(),
priority: Priority::Hard,
constraint,
}
}
#[must_use]
pub fn soft(name: impl Into<String>, weight: f32, constraint: Constraint) -> Self {
Self {
name: name.into(),
priority: Priority::Soft(weight),
constraint,
}
}
}
#[derive(Debug, Clone)]
pub struct Violation {
pub law_name: String,
pub priority: Priority,
pub residual: f32,
pub point: Vec3,
pub region: Vec3Interval,
pub evidence: Evidence,
}
#[derive(Debug, Clone, PartialEq)]
#[non_exhaustive]
pub enum UnresolvedReason {
SurfaceProximity {
radius: f32,
},
SignUndecided,
GapUnbracketed {
upper: f32,
},
GradientUnwitnessed {
claimed: f32,
worst_sampled: f32,
},
ReachabilityUndecided {
undecided_cells: usize,
},
SurfaceRatioUnbracketed {
lo: f32,
hi: f32,
},
TemperatureUnbracketed {
lo_c: f32,
hi_c: f32,
},
ThermalGridUnusable {
why: &'static str,
},
VolumeUnbracketed {
lo: f32,
hi: f32,
},
}
#[derive(Debug, Clone)]
pub struct Unresolved {
pub law_name: String,
pub priority: Priority,
pub point: Vec3,
pub region: Vec3Interval,
pub reason: UnresolvedReason,
}
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum HardVerdict {
Proven,
Violated,
Undecided,
}
impl HardVerdict {
#[must_use]
pub const fn is_proven(self) -> bool {
matches!(self, Self::Proven)
}
}
#[derive(Debug, Clone)]
pub struct LawReport {
pub total_laws: usize,
pub passed: usize,
pub violations: Vec<Violation>,
pub unresolved: Vec<Unresolved>,
}
impl LawReport {
#[must_use]
pub const fn all_passed(&self) -> bool {
self.violations.is_empty() && self.unresolved.is_empty()
}
#[must_use]
pub fn has_hard_violations(&self) -> bool {
self.violations.iter().any(|v| v.priority == Priority::Hard)
}
#[must_use]
pub fn proven_violations(&self) -> Vec<&Violation> {
self.violations
.iter()
.filter(|v| v.evidence.is_proof())
.collect()
}
#[must_use]
pub const fn has_unresolved(&self) -> bool {
!self.unresolved.is_empty()
}
#[must_use]
pub fn hard_verdict(&self) -> HardVerdict {
if self.has_hard_violations() {
HardVerdict::Violated
} else if self.unresolved.iter().any(|u| u.priority == Priority::Hard) {
HardVerdict::Undecided
} else {
HardVerdict::Proven
}
}
}
struct Verdict {
violation: Option<Violation>,
unresolved: Option<Unresolved>,
}
pub const BALL_PROBE_DEPTH: u32 = 4;
const BISECT_ITERS: u32 = 12;
#[derive(Debug, Clone)]
pub struct CheckConfig {
pub aabb_min: Vec3,
pub aabb_max: Vec3,
pub resolution: usize,
}
impl Default for CheckConfig {
fn default() -> Self {
Self {
aabb_min: Vec3::splat(-5.0),
aabb_max: Vec3::splat(5.0),
resolution: 8,
}
}
}
fn sdf_eval(node: &SdfNode, point: Vec3) -> f32 {
alice_sdf::eval(node, point)
}
fn sdf_interval(node: &SdfNode, bounds: Vec3Interval) -> Interval {
alice_sdf::interval::eval_interval(node, bounds)
}
fn cube(c: Vec3, h: f32) -> Vec3Interval {
Vec3Interval {
x: Interval::new(c.x - h, c.x + h),
y: Interval::new(c.y - h, c.y + h),
z: Interval::new(c.z - h, c.z + h),
}
}
const fn box_center(b: Vec3Interval) -> Vec3 {
Vec3::new(
b.x.lo.midpoint(b.x.hi),
b.y.lo.midpoint(b.y.hi),
b.z.lo.midpoint(b.z.hi),
)
}
fn box_expand(b: Vec3Interval, amount: f32) -> Vec3Interval {
Vec3Interval {
x: b.x.expand(amount),
y: b.y.expand(amount),
z: b.z.expand(amount),
}
}
const fn box_nearest(b: Vec3Interval, p: Vec3) -> Vec3 {
Vec3::new(
p.x.clamp(b.x.lo, b.x.hi),
p.y.clamp(b.y.lo, b.y.hi),
p.z.clamp(b.z.lo, b.z.hi),
)
}
fn box_children(b: Vec3Interval) -> [Vec3Interval; 8] {
let c = box_center(b);
let split = |iv: Interval, m: f32, hi: bool| {
if hi {
Interval::new(m, iv.hi)
} else {
Interval::new(iv.lo, m)
}
};
let mut out = [b; 8];
for (i, child) in out.iter_mut().enumerate() {
*child = Vec3Interval {
x: split(b.x, c.x, i & 1 != 0),
y: split(b.y, c.y, i & 2 != 0),
z: split(b.z, c.z, i & 4 != 0),
};
}
out
}
fn interval_sign(iv: Interval) -> Option<bool> {
if iv.lo >= 0.0 {
Some(true)
} else if iv.hi < 0.0 {
Some(false)
} else {
None
}
}
fn bisect_crossing(node: &SdfNode, p: Vec3, outside_p: bool, mut q: Vec3) -> f32 {
let mut a = p;
for _ in 0..BISECT_ITERS {
let m = (a + q) * 0.5;
if (sdf_eval(node, m) >= 0.0) == outside_p {
a = m;
} else {
q = m;
}
}
p.distance(q)
}
#[derive(Debug)]
enum BallProbe {
Clear,
Crossing { distance: f32 },
Undecided { region: Vec3Interval },
}
fn probe_ball(node: &SdfNode, centre: Vec3, radius: f32) -> BallProbe {
let f_centre = sdf_eval(node, centre);
if f_centre == 0.0 {
return BallProbe::Crossing { distance: 0.0 };
}
let outside = f_centre > 0.0;
let mut best: Option<f32> = None;
let mut undecided: Option<Vec3Interval> = None;
let l = eval_lipschitz(node);
let mut stack: Vec<(Vec3Interval, u32)> = vec![(cube(centre, radius), 0)];
while let Some((bx, depth)) = stack.pop() {
let near = box_nearest(bx, centre);
let near_dist = centre.distance(near);
if near_dist > radius {
continue; }
if best.is_some_and(|d| near_dist >= d) {
continue; }
let iv = sdf_interval(node, bx);
let sign = match refine(node, bx, iv, l) {
Refined::Enclosure(tight) => interval_sign(tight).or_else(|| interval_sign(iv)),
Refined::Contradiction => None,
};
match sign {
Some(sign) if sign == outside => {}
Some(_) => {
let d = bisect_crossing(node, centre, outside, near);
if best.is_none_or(|bd| d < bd) {
best = Some(d);
}
}
None if depth >= BALL_PROBE_DEPTH => {
let leaf_centre = box_center(bx);
if centre.distance(leaf_centre) <= radius
&& (sdf_eval(node, leaf_centre) >= 0.0) != outside
{
let d = bisect_crossing(node, centre, outside, leaf_centre);
if best.is_none_or(|bd| d < bd) {
best = Some(d);
}
} else if undecided.is_none() {
undecided = Some(bx);
}
}
None => {
let mut children = box_children(bx);
children.sort_by(|x, y| {
let dx = centre.distance(box_nearest(*x, centre));
let dy = centre.distance(box_nearest(*y, centre));
dy.partial_cmp(&dx).unwrap_or(std::cmp::Ordering::Equal)
});
for child in children {
stack.push((child, depth + 1));
}
}
}
}
match (best, undecided) {
(Some(distance), _) => BallProbe::Crossing { distance },
(None, Some(region)) => BallProbe::Undecided { region },
(None, None) => BallProbe::Clear,
}
}
fn surface_distance(node: &SdfNode, p: Vec3, floor: f32) -> Option<f32> {
let field = sdf_eval(node, p).abs();
let cap = crate::det_math::max(
crate::det_math::max(floor, field * SURFACE_SEARCH_SLACK),
f32::EPSILON,
);
let mut r = crate::det_math::max(field, cap * 1.0e-3);
for _ in 0..SURFACE_SEARCH_STEPS {
if let BallProbe::Crossing { distance } = probe_ball(node, p, crate::det_math::min(r, cap))
{
return Some(distance);
}
if r >= cap {
break;
}
r *= 2.0;
}
None
}
const SURFACE_SEARCH_STEPS: u32 = 10;
const SURFACE_SEARCH_SLACK: f32 = 16.0;
fn search_cap(config: &CheckConfig) -> f32 {
(config.aabb_max - config.aabb_min).length()
}
fn accumulate_thickness(
node: &SdfNode,
center: Vec3,
bounds: Vec3Interval,
required: f32,
worst: &mut Option<(f32, Vec3, Vec3Interval)>,
undecided: &mut Option<(Vec3, Vec3Interval, f32)>,
) {
match probe_ball(node, center, required) {
BallProbe::Clear => {}
BallProbe::Crossing { distance } => {
let residual = distance - required; match worst {
Some((w, _, _)) if residual >= *w => {}
_ => *worst = Some((residual, center, bounds)),
}
}
BallProbe::Undecided { region } => {
if undecided.is_none() {
*undecided = Some((center, region, required));
}
}
}
}
fn thickness_verdict(
law_name: &str,
priority: Priority,
evidence: Evidence,
worst: Option<(f32, Vec3, Vec3Interval)>,
undecided: Option<(Vec3, Vec3Interval, f32)>,
) -> Verdict {
Verdict {
violation: worst.map(|(residual, point, region)| Violation {
law_name: law_name.to_string(),
priority,
residual,
point,
region,
evidence,
}),
unresolved: undecided.map(|(point, region, radius)| Unresolved {
law_name: law_name.to_string(),
priority,
point,
region,
reason: UnresolvedReason::SurfaceProximity { radius },
}),
}
}
enum PairProbe {
Clear,
Witness { point: Vec3, residual: f32 },
Undecided { region: Vec3Interval },
}
fn lipschitz_enclosure(bx: Vec3Interval, f_centre: f32, l: f32) -> Option<Interval> {
if !l.is_finite() || l <= 0.0 || !f_centre.is_finite() || f_centre < 0.0 {
return None;
}
let (dx, dy, dz) = (bx.x.hi - bx.x.lo, bx.y.hi - bx.y.lo, bx.z.hi - bx.z.lo);
let rho = 0.5 * dx.mul_add(dx, dy.mul_add(dy, dz * dz)).sqrt();
if !rho.is_finite() {
return None;
}
let slack = l * rho;
Some(Interval::new(f_centre - slack, f_centre + slack))
}
#[derive(Debug, Clone, Copy)]
enum Refined {
Enclosure(Interval),
Contradiction,
}
fn refine_with_lipschitz(bx: Vec3Interval, iv: Interval, f_centre: f32, l: f32) -> Refined {
let Some(e) = lipschitz_enclosure(bx, f_centre, l) else {
return Refined::Enclosure(iv);
};
let (lo, hi) = (
crate::det_math::max(iv.lo, e.lo),
crate::det_math::min(iv.hi, e.hi),
);
if lo > hi {
return Refined::Contradiction;
}
Refined::Enclosure(Interval::new(lo, hi))
}
fn refine(node: &SdfNode, bx: Vec3Interval, iv: Interval, l: f32) -> Refined {
if !l.is_finite() || l <= 0.0 {
return Refined::Enclosure(iv);
}
refine_with_lipschitz(bx, iv, sdf_eval(node, box_center(bx)), l)
}
fn probe_pair(
a: &SdfNode,
b: &SdfNode,
cell: Vec3Interval,
decided_ok: &dyn Fn(Interval, Interval) -> bool,
witness: &dyn Fn(f32, f32) -> bool,
residual: &dyn Fn(f32, f32) -> f32,
) -> PairProbe {
let mut worst: Option<(f32, Vec3)> = None;
let mut undecided: Option<Vec3Interval> = None;
let mut stack: Vec<(Vec3Interval, u32)> = vec![(cell, 0)];
let (l_a, l_b) = (eval_lipschitz(a), eval_lipschitz(b));
while let Some((bx, depth)) = stack.pop() {
let ia = sdf_interval(a, bx);
let ib = sdf_interval(b, bx);
let c = box_center(bx);
let (fa, fb) = (sdf_eval(a, c), sdf_eval(b, c));
let (Refined::Enclosure(ta), Refined::Enclosure(tb)) = (
refine_with_lipschitz(bx, ia, fa, l_a),
refine_with_lipschitz(bx, ib, fb, l_b),
) else {
if undecided.is_none() {
undecided = Some(bx);
}
continue;
};
if decided_ok(ia, ib) {
continue;
}
if witness(fa, fb) {
let r = raw_min_f32(
residual(fa, fb),
-f32::EPSILON,
crate::raw_minmax::Reason::NonzeroConst,
);
if worst.is_none_or(|(w, _)| r < w) {
worst = Some((r, c));
}
}
if decided_ok(ta, tb) {
continue;
}
if depth >= BALL_PROBE_DEPTH {
if undecided.is_none() {
undecided = Some(bx);
}
} else {
for child in box_children(bx) {
stack.push((child, depth + 1));
}
}
}
match (worst, undecided) {
(Some((residual, point)), _) => PairProbe::Witness { point, residual },
(None, Some(region)) => PairProbe::Undecided { region },
(None, None) => PairProbe::Clear,
}
}
fn pair_verdict(
law_name: &str,
priority: Priority,
worst: Option<(f32, Vec3, Vec3Interval)>,
undecided: Option<(Vec3, Vec3Interval)>,
) -> Verdict {
Verdict {
violation: worst.map(|(residual, point, region)| Violation {
law_name: law_name.to_string(),
priority,
residual,
point,
region,
evidence: Evidence::Witnessed,
}),
unresolved: undecided.map(|(point, region)| Unresolved {
law_name: law_name.to_string(),
priority,
point,
region,
reason: UnresolvedReason::SignUndecided,
}),
}
}
fn check_non_overlap(
a: &SdfNode,
b: &SdfNode,
law_name: &str,
priority: Priority,
config: &CheckConfig,
) -> Verdict {
let mut worst: Option<(f32, Vec3, Vec3Interval)> = None;
let mut undecided: Option<(Vec3, Vec3Interval)> = None;
for (center, bounds) in GridSampler::new(config) {
match probe_pair(
a,
b,
bounds,
&|ia, ib| ia.lo >= 0.0 || ib.lo >= 0.0,
&|fa, fb| fa < 0.0 && fb < 0.0,
&|fa, fb| crate::det_math::max(fa, fb),
) {
PairProbe::Clear => {}
PairProbe::Witness { point, residual } => match &worst {
Some((w, _, _)) if residual >= *w => {}
_ => worst = Some((residual, point, bounds)),
},
PairProbe::Undecided { region } => {
if undecided.is_none() {
undecided = Some((center, region));
}
}
}
}
let cap = search_cap(config);
let worst = worst.map(|(field, point, region)| {
let depth = surface_distance(a, point, cap)
.into_iter()
.chain(surface_distance(b, point, cap))
.fold(f32::INFINITY, crate::det_math::min);
let residual = if depth.is_finite() {
-depth
} else {
field };
(
raw_min_f32(
residual,
-f32::EPSILON,
crate::raw_minmax::Reason::NonzeroConst,
),
point,
region,
)
});
pair_verdict(law_name, priority, worst, undecided)
}
fn check_containment(
inner: &SdfNode,
outer: &SdfNode,
law_name: &str,
priority: Priority,
config: &CheckConfig,
) -> Verdict {
let mut worst: Option<(f32, Vec3, Vec3Interval)> = None;
let mut undecided: Option<(Vec3, Vec3Interval)> = None;
for (center, bounds) in GridSampler::new(config) {
match probe_pair(
inner,
outer,
bounds,
&|ii, io| ii.lo >= 0.0 || io.hi <= 0.0,
&|fi, fo| fi < 0.0 && fo > 0.0,
&|_, fo| -fo,
) {
PairProbe::Clear => {}
PairProbe::Witness { point, residual } => match &worst {
Some((w, _, _)) if residual >= *w => {}
_ => worst = Some((residual, point, bounds)),
},
PairProbe::Undecided { region } => {
if undecided.is_none() {
undecided = Some((center, region));
}
}
}
}
let cap = search_cap(config);
let worst = worst.map(|(field, point, region)| {
let residual = surface_distance(outer, point, cap).map_or(field, |d| -d);
(
raw_min_f32(
residual,
-f32::EPSILON,
crate::raw_minmax::Reason::NonzeroConst,
),
point,
region,
)
});
pair_verdict(law_name, priority, worst, undecided)
}
fn check_min_thickness(
node: &SdfNode,
min_thickness: f32,
law_name: &str,
priority: Priority,
config: &CheckConfig,
) -> Verdict {
let mut worst = None;
let mut undecided = None;
for (center, bounds) in GridSampler::new(config) {
if sdf_eval(node, center) >= 0.0 {
continue; }
accumulate_thickness(
node,
center,
bounds,
min_thickness,
&mut worst,
&mut undecided,
);
}
thickness_verdict(law_name, priority, Evidence::Witnessed, worst, undecided)
}
fn check_stress(
node: &SdfNode,
load_points: &[(Vec3, f32)],
min_thickness_factor: f32,
law_name: &str,
priority: Priority,
config: &CheckConfig,
) -> Verdict {
let mut worst = None;
let mut undecided = None;
for (center, bounds) in GridSampler::new(config) {
if sdf_eval(node, center) >= 0.0 {
continue; }
for &(lp, force) in load_points {
let search_radius = raw_max_f32(force, 1.0, crate::raw_minmax::Reason::NonzeroConst);
if center.distance(lp) > search_radius {
continue;
}
let required = force * min_thickness_factor;
if required <= 0.0 {
continue;
}
accumulate_thickness(node, center, bounds, required, &mut worst, &mut undecided);
}
}
thickness_verdict(
law_name,
priority,
Evidence::Modelled {
model: MODEL_STRESS,
},
worst,
undecided,
)
}
fn gap_exceeds(a: &SdfNode, b: &SdfNode, m: f32, config: &CheckConfig) -> bool {
let half = m * 0.5;
let (l_a, l_b) = (eval_lipschitz(a), eval_lipschitz(b));
let provably_outside = |node: &SdfNode, e: Vec3Interval, l: f32| {
let iv = sdf_interval(node, e);
matches!(refine(node, e, iv, l), Refined::Enclosure(tight) if tight.lo > 0.0)
};
for (_, cell) in GridSampler::new(config) {
let mut stack: Vec<(Vec3Interval, u32)> = vec![(cell, 0)];
while let Some((bx, depth)) = stack.pop() {
let e = box_expand(bx, half);
if provably_outside(a, e, l_a) || provably_outside(b, e, l_b) {
continue;
}
if depth >= BALL_PROBE_DEPTH {
return false;
}
for child in box_children(bx) {
stack.push((child, depth + 1));
}
}
}
true
}
fn contact_upper_bound(
a: &SdfNode,
b: &SdfNode,
min_distance: f32,
max_distance: f32,
config: &CheckConfig,
) -> Option<(f32, Vec3, Vec3Interval)> {
#[allow(clippy::cast_precision_loss)]
let cell_diag =
((config.aabb_max - config.aabb_min) / config.resolution.max(1) as f32).length();
let radius = crate::det_math::max(max_distance, min_distance) + cell_diag;
let mut upper: Option<(f32, Vec3, Vec3Interval)> = None;
for (center, bounds) in GridSampler::new(config) {
let (fa, fb) = (sdf_eval(a, center), sdf_eval(b, center));
if fa <= 0.0 || fb <= 0.0 {
continue;
}
let e = box_expand(bounds, radius);
if sdf_interval(a, e).lo > 0.0 || sdf_interval(b, e).lo > 0.0 {
continue;
}
let (BallProbe::Crossing { distance: da }, BallProbe::Crossing { distance: db }) =
(probe_ball(a, center, radius), probe_ball(b, center, radius))
else {
continue;
};
let ub = da + db;
match &upper {
Some((u, _, _)) if ub >= *u => {}
_ => upper = Some((ub, center, bounds)),
}
}
upper
}
fn check_contact(
a: &SdfNode,
b: &SdfNode,
min_distance: f32,
max_distance: f32,
law_name: &str,
priority: Priority,
config: &CheckConfig,
) -> Verdict {
let interference = check_non_overlap(a, b, law_name, priority, config);
if interference.violation.is_some() {
return interference;
}
let upper = contact_upper_bound(a, b, min_distance, max_distance, config);
let make_violation =
|residual: f32, point: Vec3, region: Vec3Interval, evidence: Evidence| Violation {
law_name: law_name.to_string(),
priority,
residual,
point,
region,
evidence,
};
if let Some((ub, point, region)) = upper {
if ub < min_distance {
return Verdict {
violation: Some(make_violation(
ub - min_distance,
point,
region,
Evidence::Witnessed,
)),
unresolved: None,
};
}
}
if gap_exceeds(a, b, max_distance, config) {
let (residual, point, region) = upper.map_or_else(
|| {
(
f32::NEG_INFINITY,
config.aabb_min,
Vec3Interval::from_bounds(config.aabb_min, config.aabb_max),
)
},
|(ub, p, r)| (max_distance - ub, p, r),
);
return Verdict {
violation: Some(make_violation(residual, point, region, Evidence::Proved)),
unresolved: None,
};
}
if let Some((ub, _, _)) = upper {
if ub <= max_distance && gap_exceeds(a, b, min_distance, config) {
return interference; }
}
Verdict {
violation: None,
unresolved: Some(Unresolved {
law_name: law_name.to_string(),
priority,
point: upper.map_or(config.aabb_min, |(_, p, _)| p),
region: upper.map_or_else(|| cube(config.aabb_min, 0.0), |(_, _, r)| r),
reason: UnresolvedReason::GapUnbracketed {
upper: upper.map_or(f32::INFINITY, |(u, _, _)| u),
},
}),
}
}
struct GridSampler {
aabb_min: Vec3,
step: Vec3,
half_step: Vec3,
n: usize,
ix: usize,
iy: usize,
iz: usize,
}
impl GridSampler {
fn new(config: &CheckConfig) -> Self {
let n = config.resolution;
let extent = config.aabb_max - config.aabb_min;
#[allow(clippy::cast_precision_loss)]
let step = extent / n as f32;
Self {
aabb_min: config.aabb_min,
step,
half_step: step * 0.5,
n,
ix: 0,
iy: 0,
iz: 0,
}
}
}
impl Iterator for GridSampler {
type Item = (Vec3, Vec3Interval);
fn next(&mut self) -> Option<Self::Item> {
if self.iz >= self.n {
return None;
}
#[allow(clippy::cast_precision_loss)]
let lo =
self.aabb_min + self.step * Vec3::new(self.ix as f32, self.iy as f32, self.iz as f32);
let center = lo + self.half_step;
let hi = lo + self.step;
let bounds = Vec3Interval {
x: Interval { lo: lo.x, hi: hi.x },
y: Interval { lo: lo.y, hi: hi.y },
z: Interval { lo: lo.z, hi: hi.z },
};
self.ix += 1;
if self.ix >= self.n {
self.ix = 0;
self.iy += 1;
if self.iy >= self.n {
self.iy = 0;
self.iz += 1;
}
}
Some((center, bounds))
}
}
fn check_one(law: &Law, config: &CheckConfig) -> Verdict {
match &law.constraint {
Constraint::NonOverlap { a, b } => check_non_overlap(a, b, &law.name, law.priority, config),
Constraint::Containment { inner, outer } => {
check_containment(inner, outer, &law.name, law.priority, config)
}
Constraint::MinThickness {
node,
min_thickness,
} => check_min_thickness(node, *min_thickness, &law.name, law.priority, config),
Constraint::Stress {
node,
load_points,
min_thickness_factor,
} => check_stress(
node,
load_points,
*min_thickness_factor,
&law.name,
law.priority,
config,
),
Constraint::Thermal {
node,
heat_sources,
search_radius,
min_surface_ratio,
} => check_thermal(
node,
heat_sources,
*search_radius,
*min_surface_ratio,
&law.name,
law.priority,
config,
),
Constraint::Contact {
a,
b,
min_distance,
max_distance,
} => check_contact(
a,
b,
*min_distance,
*max_distance,
&law.name,
law.priority,
config,
),
Constraint::Continuity { node, seed_point } => {
check_continuity(node, *seed_point, &law.name, law.priority, config)
}
Constraint::GradientBound {
node,
max_gradient,
probe,
} => check_gradient_bound(node, *max_gradient, *probe, &law.name, law.priority, config),
Constraint::Reachable { node, from, to } => {
check_reachable(node, *from, *to, &law.name, law.priority, config)
}
Constraint::VolumeConservation {
before,
after,
relative_tolerance,
} => check_volume_conservation(
before,
after,
*relative_tolerance,
&law.name,
law.priority,
config,
),
#[cfg(feature = "physics")]
Constraint::ThermalField {
node,
sources,
material,
metres_per_unit,
ambient_c,
max_temperature_c,
duration_s,
} => check_thermal_field(
node,
sources,
material,
*metres_per_unit,
*ambient_c,
*max_temperature_c,
*duration_s,
&law.name,
law.priority,
config,
),
}
}
#[must_use]
pub fn check_laws(laws: &[Law], config: &CheckConfig) -> LawReport {
let mut violations = Vec::new();
let mut unresolved = Vec::new();
for law in laws {
let verdict = check_one(law, config);
if let Some(v) = verdict.violation {
violations.push(v);
}
if let Some(u) = verdict.unresolved {
unresolved.push(u);
}
}
violations.sort_by(|a, b| {
b.residual
.abs()
.partial_cmp(&a.residual.abs())
.unwrap_or(std::cmp::Ordering::Equal)
});
let mut undecided_names: Vec<&str> = violations.iter().map(|v| v.law_name.as_str()).collect();
undecided_names.extend(unresolved.iter().map(|u| u.law_name.as_str()));
undecided_names.sort_unstable();
undecided_names.dedup();
let passed = laws.len().saturating_sub(undecided_names.len());
LawReport {
total_laws: laws.len(),
passed,
violations,
unresolved,
}
}
fn surface_ratio_bracket(
node: &SdfNode,
source: Vec3,
search_radius: f32,
config: &CheckConfig,
signs: &[CellSign],
) -> Option<(f32, f32, (Vec3, Vec3Interval))> {
let n = config.resolution;
let step = grid_step(config);
let (mut surface_lo, mut surface_hi) = (0usize, 0usize);
let (mut interior_lo, mut interior_hi) = (0usize, 0usize);
let mut representative: Option<(Vec3, Vec3Interval)> = None;
for (ix, iy, iz, idx) in grid_indices(n) {
let center = cell_center(config, step, ix, iy, iz);
if center.distance(source) > search_radius {
continue;
}
let region = cell_box(config, step, ix, iy, iz);
match signs[idx] {
CellSign::Inside => {
interior_lo += 1;
interior_hi += 1;
if representative.is_none() {
representative = Some((center, region));
}
}
CellSign::Outside => {}
CellSign::Undecided => {
surface_hi += 1;
interior_hi += 1;
if cell_straddles_surface(node, region) {
surface_lo += 1;
}
}
}
}
if interior_hi == 0 {
return None; }
let representative = representative.or_else(|| {
Some((source, cube(source, step.length() * 0.5)))
})?;
#[allow(clippy::cast_precision_loss)]
let lo = surface_lo as f32 / interior_hi as f32;
#[allow(clippy::cast_precision_loss)]
let hi = if interior_lo == 0 {
f32::INFINITY } else {
surface_hi as f32 / interior_lo as f32
};
Some((lo, hi, representative))
}
fn cell_straddles_surface(node: &SdfNode, region: Vec3Interval) -> bool {
let lo = Vec3::new(region.x.lo, region.y.lo, region.z.lo);
let hi = Vec3::new(region.x.hi, region.y.hi, region.z.hi);
let first = sdf_eval(node, lo) >= 0.0;
for c in 1..8u32 {
let p = Vec3::new(
if c & 1 == 0 { lo.x } else { hi.x },
if c & 2 == 0 { lo.y } else { hi.y },
if c & 4 == 0 { lo.z } else { hi.z },
);
if (sdf_eval(node, p) >= 0.0) != first {
return true;
}
}
(sdf_eval(node, box_center(region)) >= 0.0) != first
}
fn check_thermal(
node: &SdfNode,
heat_sources: &[Vec3],
search_radius: f32,
min_surface_ratio: f32,
law_name: &str,
priority: Priority,
config: &CheckConfig,
) -> Verdict {
let step = grid_step(config);
let signs = classify_cells(node, config, step);
let mut worst: Option<(f32, Vec3, Vec3Interval)> = None;
let mut undecided: Option<(Vec3, Vec3Interval, f32, f32)> = None;
for &source in heat_sources {
let Some((lo, hi, (point, region))) =
surface_ratio_bracket(node, source, search_radius, config, &signs)
else {
continue;
};
if hi < min_surface_ratio {
let residual = hi - min_surface_ratio; match &worst {
Some((w, _, _)) if residual >= *w => {}
_ => worst = Some((residual, point, region)),
}
} else if lo < min_surface_ratio && undecided.is_none() {
undecided = Some((point, region, lo, hi));
}
}
Verdict {
violation: worst.map(|(residual, point, region)| Violation {
law_name: law_name.to_string(),
priority,
residual,
point,
region,
evidence: Evidence::Modelled {
model: MODEL_THERMAL,
},
}),
unresolved: undecided.map(|(point, region, lo, hi)| Unresolved {
law_name: law_name.to_string(),
priority,
point,
region,
reason: UnresolvedReason::SurfaceRatioUnbracketed { lo, hi },
}),
}
}
#[allow(dead_code)]
const FRACTION_SAMPLES: usize = 512;
#[allow(dead_code)]
fn material_fraction(node: &SdfNode, p: Vec3, r: f32) -> f32 {
const GOLDEN_ANGLE: f32 = 2.399_963_2; let mut inside = 0_u32;
for i in 0..FRACTION_SAMPLES {
#[allow(clippy::cast_precision_loss)]
let t = (i as f32 + 0.5) / FRACTION_SAMPLES as f32;
let z = 2.0f32.mul_add(-t, 1.0);
let rho = crate::det_math::max(1.0 - z * z, 0.0).sqrt();
#[allow(clippy::cast_precision_loss)]
let theta = GOLDEN_ANGLE * i as f32;
let (sin_t, cos_t) = alice_det_math::sin_cos(theta);
let dir = Vec3::new(rho * cos_t, rho * sin_t, z);
if sdf_eval(node, p + dir * r) < 0.0 {
inside += 1;
}
}
#[allow(clippy::cast_precision_loss)]
let frac = f32::from(u16::try_from(inside).unwrap_or(u16::MAX)) / FRACTION_SAMPLES as f32;
frac
}
#[cfg(feature = "physics")]
const THERMAL_MAX_STEPS: u32 = 6000;
#[cfg(feature = "physics")]
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
enum SurfaceCooling {
Adiabatic,
Isothermal,
}
#[cfg(feature = "physics")]
struct ThermalScene<'a> {
solid: &'a [bool],
sources: &'a [(Vec3, f32)],
material: &'a alice_physics::transient_thermal::ThermalMaterial,
metres_per_unit: f32,
ambient_c: f32,
duration_s: f32,
}
#[cfg(feature = "physics")]
fn solve_peak_temperature(
scene: &ThermalScene<'_>,
cooling: SurfaceCooling,
config: &CheckConfig,
) -> Option<f32> {
use alice_physics::transient_thermal::{stable_dt_3d, transient_step_3d};
let ThermalScene {
solid,
sources,
material,
metres_per_unit,
ambient_c,
duration_s,
} = *scene;
let n = config.resolution;
let step = grid_step(config);
let dx_m = step.x * metres_per_unit;
let cell_volume = dx_m * dx_m * dx_m;
let mut t = vec![ambient_c; n * n * n];
let mut source_cells: Vec<(usize, f32)> = Vec::new();
for &(p, watts) in sources {
let rel = (p - config.aabb_min) / step;
if rel.x < 0.0 || rel.y < 0.0 || rel.z < 0.0 {
continue;
}
let (ix, iy, iz) = (
grid_index(rel.x, n),
grid_index(rel.y, n),
grid_index(rel.z, n),
);
let idx = ix + iy * n + iz * n * n;
if solid[idx] {
source_cells.push((idx, watts));
}
}
let dt_max = stable_dt_3d(&t, material, dx_m);
if !dt_max.is_finite() || dt_max <= 0.0 {
return None;
}
let dt = crate::det_math::min(dt_max * 0.5, duration_s);
#[allow(clippy::cast_possible_truncation, clippy::cast_sign_loss)]
let steps = (duration_s / dt).ceil() as u32;
if steps > THERMAL_MAX_STEPS {
return None;
}
let neighbours = |idx: usize| -> [Option<usize>; 6] {
let (ix, iy, iz) = (idx % n, (idx / n) % n, idx / (n * n));
[
ix.checked_sub(1).map(|v| v + iy * n + iz * n * n),
(ix + 1 < n).then_some((ix + 1) + iy * n + iz * n * n),
iy.checked_sub(1).map(|v| ix + v * n + iz * n * n),
(iy + 1 < n).then_some(ix + (iy + 1) * n + iz * n * n),
iz.checked_sub(1).map(|v| ix + iy * n + v * n * n),
(iz + 1 < n).then_some(ix + iy * n + (iz + 1) * n * n),
]
};
for _ in 0..steps {
for &(idx, watts) in &source_cells {
let rho_cp = material.heat_capacity_at(t[idx]);
if rho_cp > 0.0 && rho_cp.is_finite() {
t[idx] += watts * dt / (rho_cp * cell_volume);
}
}
transient_step_3d(&mut t, n, n, n, material, dx_m, dt);
match cooling {
SurfaceCooling::Isothermal => {
for (idx, &is_solid) in solid.iter().enumerate() {
if !is_solid {
t[idx] = ambient_c;
}
}
}
SurfaceCooling::Adiabatic => {
let snapshot = t.clone();
for idx in 0..solid.len() {
if solid[idx] {
continue;
}
let mut sum = 0.0_f32;
let mut count = 0_u32;
for nb in neighbours(idx).into_iter().flatten() {
if solid[nb] {
sum += snapshot[nb];
count += 1;
}
}
t[idx] = if count == 0 {
ambient_c
} else {
#[allow(clippy::cast_precision_loss)]
let c = f32::from(u16::try_from(count).unwrap_or(u16::MAX));
sum / c
};
}
}
}
}
let peak = t
.iter()
.zip(solid)
.filter_map(|(v, &is_solid)| is_solid.then_some(*v))
.fold(f32::NEG_INFINITY, crate::det_math::max);
peak.is_finite().then_some(peak)
}
#[cfg(feature = "physics")]
fn thermal_source_placement_problem(
sources: &[(Vec3, f32)],
solid: &[bool],
loose: &[bool],
config: &CheckConfig,
step: Vec3,
) -> Option<&'static str> {
let n = config.resolution;
for &(p, watts) in sources {
if watts == 0.0 {
continue; }
let rel = (p - config.aabb_min) / step;
if rel.x < 0.0 || rel.y < 0.0 || rel.z < 0.0 {
return Some("熱源が検査 AABB の外にある");
}
let (ix, iy, iz) = (
grid_index(rel.x, n),
grid_index(rel.y, n),
grid_index(rel.z, n),
);
let start = (ix, iy, iz, ix + iy * n + iz * n * n);
if !solid[start.3] {
return Some("熱源が内部と確定したセルに乗っていない (解像度を上げる)");
}
let strict = flood_mask(n, start, &|i| solid[i]);
let relaxed = flood_mask(n, start, &|i| loose[i]);
if (0..solid.len()).any(|i| solid[i] && relaxed[i] && !strict[i]) {
return Some(
"熱源のある材料が格子の都合で分断されている \
(形状としては繋がっているので解像度を上げる)",
);
}
}
None
}
#[cfg(feature = "physics")]
#[allow(clippy::too_many_arguments)]
fn check_thermal_field(
node: &SdfNode,
sources: &[(Vec3, f32)],
material: &alice_physics::transient_thermal::ThermalMaterial,
metres_per_unit: f32,
ambient_c: f32,
max_temperature_c: f32,
duration_s: f32,
law_name: &str,
priority: Priority,
config: &CheckConfig,
) -> Verdict {
let step = grid_step(config);
let unusable = |why: &'static str| Verdict {
violation: None,
unresolved: Some(Unresolved {
law_name: law_name.to_string(),
priority,
point: box_center(Vec3Interval::from_bounds(config.aabb_min, config.aabb_max)),
region: Vec3Interval::from_bounds(config.aabb_min, config.aabb_max),
reason: UnresolvedReason::ThermalGridUnusable { why },
}),
};
if config.resolution < 3 {
return unusable("格子が各軸 3 セル未満 (transient_step_3d が解かない)");
}
let aniso = raw_max_f32(
(step.x - step.y).abs(),
(step.x - step.z).abs(),
crate::raw_minmax::Reason::BothAbs,
);
if aniso > step.x.abs() * 1.0e-4 {
return unusable("セルが立方でない (検査 AABB を立方体にする)");
}
if !(metres_per_unit > 0.0 && duration_s > 0.0) {
return unusable("metres_per_unit / duration_s が正でない");
}
let signs = classify_cells(node, config, step);
let solid: Vec<bool> = signs.iter().map(|s| *s == CellSign::Inside).collect();
let loose: Vec<bool> = signs.iter().map(|s| *s != CellSign::Outside).collect();
if !solid.iter().any(|b| *b) {
return unusable("内部と確定したセルが 1 つも無い (解像度を上げる)");
}
if let Some(why) = thermal_source_placement_problem(sources, &solid, &loose, config, step) {
return unusable(why);
}
let scene = ThermalScene {
solid: &solid,
sources,
material,
metres_per_unit,
ambient_c,
duration_s,
};
let solve = |cooling| solve_peak_temperature(&scene, cooling, config);
let (Some(hi_c), Some(lo_c)) = (
solve(SurfaceCooling::Adiabatic),
solve(SurfaceCooling::Isothermal),
) else {
return unusable("CFL 条件で step 数が上限を超えた (解像度を下げるか時間を短く)");
};
if lo_c > max_temperature_c {
return Verdict {
violation: Some(Violation {
law_name: law_name.to_string(),
priority,
residual: max_temperature_c - lo_c, point: box_center(Vec3Interval::from_bounds(config.aabb_min, config.aabb_max)),
region: Vec3Interval::from_bounds(config.aabb_min, config.aabb_max),
evidence: Evidence::Modelled {
model: MODEL_THERMAL_FIELD,
},
}),
unresolved: None,
};
}
if hi_c <= max_temperature_c {
return Verdict {
violation: None,
unresolved: None,
};
}
Verdict {
violation: None,
unresolved: Some(Unresolved {
law_name: law_name.to_string(),
priority,
point: box_center(Vec3Interval::from_bounds(config.aabb_min, config.aabb_max)),
region: Vec3Interval::from_bounds(config.aabb_min, config.aabb_max),
reason: UnresolvedReason::TemperatureUnbracketed { lo_c, hi_c },
}),
}
}
#[allow(clippy::cast_possible_truncation, clippy::cast_sign_loss)]
fn grid_index(coord: f32, n: usize) -> usize {
(raw_max_f32(coord, 0.0, crate::raw_minmax::Reason::IntCast) as usize).min(n - 1)
}
fn flood_mask(
n: usize,
start: (usize, usize, usize, usize),
allowed: &dyn Fn(usize) -> bool,
) -> Vec<bool> {
let mut visited = vec![false; n * n * n];
if !allowed(start.3) {
return visited;
}
let mut queue = std::collections::VecDeque::new();
visited[start.3] = true;
queue.push_back((start.0, start.1, start.2));
while let Some((x, y, z)) = queue.pop_front() {
let neighbors = [
x.checked_sub(1).map(|v| (v, y, z)),
(x + 1 < n).then_some((x + 1, y, z)),
y.checked_sub(1).map(|v| (x, v, z)),
(y + 1 < n).then_some((x, y + 1, z)),
z.checked_sub(1).map(|v| (x, y, v)),
(z + 1 < n).then_some((x, y, z + 1)),
];
for (nx, ny, nz) in neighbors.into_iter().flatten() {
let nidx = nx + ny * n + nz * n * n;
if !visited[nidx] && allowed(nidx) {
visited[nidx] = true;
queue.push_back((nx, ny, nz));
}
}
}
visited
}
fn check_continuity(
node: &SdfNode,
seed_point: Vec3,
law_name: &str,
priority: Priority,
config: &CheckConfig,
) -> Verdict {
let violation = |residual: f32, point: Vec3, region: Vec3Interval| Violation {
law_name: law_name.to_string(),
priority,
residual,
point,
region,
evidence: Evidence::Modelled {
model: MODEL_CONTINUITY,
},
};
let undecided_verdict = |point: Vec3, region: Vec3Interval, cells: usize| Verdict {
violation: None,
unresolved: Some(Unresolved {
law_name: law_name.to_string(),
priority,
point,
region,
reason: UnresolvedReason::ReachabilityUndecided {
undecided_cells: cells,
},
}),
};
let seed_dist = sdf_eval(node, seed_point);
if seed_dist >= 0.0 {
return Verdict {
violation: Some(violation(seed_dist, seed_point, point_box(seed_point))),
unresolved: None,
};
}
let n = config.resolution;
let step = grid_step(config);
let signs = classify_cells(node, config, step);
let undecided_cells = signs.iter().filter(|s| **s == CellSign::Undecided).count();
let rel = (seed_point - config.aabb_min) / step;
let (sx, sy, sz) = (
grid_index(rel.x, n),
grid_index(rel.y, n),
grid_index(rel.z, n),
);
let seed = (sx, sy, sz, sx + sy * n + sz * n * n);
let seed_region = cell_box(config, step, sx, sy, sz);
let inside_total = signs.iter().filter(|s| **s == CellSign::Inside).count();
if inside_total == 0 {
return undecided_verdict(seed_point, seed_region, undecided_cells);
}
let loose = flood_mask(n, seed, &|i| signs[i] != CellSign::Outside);
let stranded =
grid_indices(n).find(|&(_, _, _, idx)| signs[idx] == CellSign::Inside && !loose[idx]);
if let Some((ix, iy, iz, _)) = stranded {
let unreachable = grid_indices(n)
.filter(|&(_, _, _, idx)| signs[idx] == CellSign::Inside && !loose[idx])
.count();
#[allow(clippy::cast_precision_loss)]
let ratio = unreachable as f32 / inside_total as f32;
let region = cell_box(config, step, ix, iy, iz);
return Verdict {
violation: Some(violation(-ratio, box_center(region), region)),
unresolved: None,
};
}
let strict = flood_mask(n, seed, &|i| signs[i] == CellSign::Inside);
let all_inside_linked = (0..signs.len()).all(|i| signs[i] != CellSign::Inside || strict[i]);
let no_detached_undecided =
(0..signs.len()).all(|i| signs[i] != CellSign::Undecided || loose[i]);
if all_inside_linked && no_detached_undecided {
return Verdict {
violation: None,
unresolved: None,
};
}
undecided_verdict(box_center(seed_region), seed_region, undecided_cells)
}
const fn point_box(p: Vec3) -> Vec3Interval {
Vec3Interval {
x: Interval { lo: p.x, hi: p.x },
y: Interval { lo: p.y, hi: p.y },
z: Interval { lo: p.z, hi: p.z },
}
}
fn cell_box(config: &CheckConfig, step: Vec3, ix: usize, iy: usize, iz: usize) -> Vec3Interval {
#[allow(clippy::cast_precision_loss)]
let lo = config.aabb_min + step * Vec3::new(ix as f32, iy as f32, iz as f32);
let hi = lo + step;
Vec3Interval {
x: Interval { lo: lo.x, hi: hi.x },
y: Interval { lo: lo.y, hi: hi.y },
z: Interval { lo: lo.z, hi: hi.z },
}
}
fn check_gradient_bound(
node: &SdfNode,
max_gradient: f32,
probe: f32,
law_name: &str,
priority: Priority,
config: &CheckConfig,
) -> Verdict {
let claimed = alice_sdf::interval::eval_lipschitz(node);
if claimed <= max_gradient {
return Verdict {
violation: None,
unresolved: None,
};
}
let n = config.resolution;
let extent = config.aabb_max - config.aabb_min;
#[allow(clippy::cast_precision_loss)]
let step = extent / (n as f32);
let h = if probe > 0.0 { probe } else { step.x * 0.25 };
let mut worst = 0.0_f32;
let mut worst_at = config.aabb_min;
let mut worst_cell = (0usize, 0usize, 0usize);
for (ix, iy, iz, _) in grid_indices(n) {
#[allow(clippy::cast_precision_loss)]
let center =
config.aabb_min + step * Vec3::new(ix as f32, iy as f32, iz as f32) + step * 0.5;
let fp = sdf_eval(node, center);
for axis in [Vec3::X, Vec3::Y, Vec3::Z] {
let q = center + axis * h;
let fq = sdf_eval(node, q);
if fp < 0.0 && fq < 0.0 {
continue;
}
if !fp.is_finite() || !fq.is_finite() {
continue;
}
let quotient = (fp - fq).abs() / h;
if quotient > worst {
worst = quotient;
worst_at = center;
worst_cell = (ix, iy, iz);
}
}
}
if worst > max_gradient {
return Verdict {
violation: Some(Violation {
law_name: law_name.to_string(),
priority,
residual: worst - max_gradient,
point: worst_at,
region: cell_box(config, step, worst_cell.0, worst_cell.1, worst_cell.2),
evidence: Evidence::Witnessed,
}),
unresolved: None,
};
}
Verdict {
violation: None,
unresolved: Some(Unresolved {
law_name: law_name.to_string(),
priority,
point: worst_at,
region: cell_box(config, step, worst_cell.0, worst_cell.1, worst_cell.2),
reason: UnresolvedReason::GradientUnwitnessed {
claimed,
worst_sampled: worst,
},
}),
}
}
fn flood_reaches(
n: usize,
a: (usize, usize, usize, usize),
b: (usize, usize, usize),
allowed: &dyn Fn(usize) -> bool,
) -> bool {
if !allowed(a.3) {
return false;
}
let mut visited = vec![false; n * n * n];
let mut queue = std::collections::VecDeque::new();
visited[a.3] = true;
queue.push_back((a.0, a.1, a.2));
while let Some((x, y, z)) = queue.pop_front() {
if (x, y, z) == b {
return true;
}
let neighbors = [
x.checked_sub(1).map(|v| (v, y, z)),
(x + 1 < n).then_some((x + 1, y, z)),
y.checked_sub(1).map(|v| (x, v, z)),
(y + 1 < n).then_some((x, y + 1, z)),
z.checked_sub(1).map(|v| (x, y, v)),
(z + 1 < n).then_some((x, y, z + 1)),
];
for (nx, ny, nz) in neighbors.into_iter().flatten() {
let nidx = nx + ny * n + nz * n * n;
if !visited[nidx] && allowed(nidx) {
visited[nidx] = true;
queue.push_back((nx, ny, nz));
}
}
}
false
}
const INTERVAL_SLACK: f32 = 1.0e-3;
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
enum CellSign {
Inside,
Outside,
Undecided,
}
fn classify_cells(node: &SdfNode, config: &CheckConfig, step: Vec3) -> Vec<CellSign> {
let n = config.resolution;
let slack =
INTERVAL_SLACK * raw_max_f32(step.length(), 1.0, crate::raw_minmax::Reason::NonzeroConst);
let mut out = vec![CellSign::Undecided; n * n * n];
for (ix, iy, iz, idx) in grid_indices(n) {
let iv = sdf_interval(node, cell_box(config, step, ix, iy, iz));
out[idx] = if iv.lo > slack {
CellSign::Outside
} else if iv.hi < -slack {
CellSign::Inside
} else {
CellSign::Undecided
};
}
out
}
fn grid_step(config: &CheckConfig) -> Vec3 {
#[allow(clippy::cast_precision_loss)]
let n = config.resolution as f32;
(config.aabb_max - config.aabb_min) / n
}
fn cell_center(config: &CheckConfig, step: Vec3, ix: usize, iy: usize, iz: usize) -> Vec3 {
#[allow(clippy::cast_precision_loss)]
let lo = config.aabb_min + step * Vec3::new(ix as f32, iy as f32, iz as f32);
lo + step * 0.5
}
fn check_reachable(
node: &SdfNode,
from: Vec3,
to: Vec3,
law_name: &str,
priority: Priority,
config: &CheckConfig,
) -> Verdict {
let n = config.resolution;
let extent = config.aabb_max - config.aabb_min;
#[allow(clippy::cast_precision_loss)]
let step = extent / (n as f32);
let slack =
INTERVAL_SLACK * raw_max_f32(step.length(), 1.0, crate::raw_minmax::Reason::NonzeroConst);
let mut proven_inside = vec![false; n * n * n];
let mut proven_outside = vec![false; n * n * n];
let mut undecided: usize = 0;
for (ix, iy, iz, idx) in grid_indices(n) {
let iv = alice_sdf::interval::eval_interval(node, cell_box(config, step, ix, iy, iz));
if iv.lo > slack {
proven_outside[idx] = true;
} else if iv.hi < -slack {
proven_inside[idx] = true;
} else {
undecided += 1;
}
}
let cell_of = |p: Vec3| -> Option<(usize, usize, usize, usize)> {
let rel = (p - config.aabb_min) / step;
if rel.x < 0.0 || rel.y < 0.0 || rel.z < 0.0 {
return None;
}
let (x, y, z) = (
grid_index(rel.x, n),
grid_index(rel.y, n),
grid_index(rel.z, n),
);
Some((x, y, z, x + y * n + z * n * n))
};
let (Some(a), Some(b)) = (cell_of(from), cell_of(to)) else {
return Verdict {
violation: None,
unresolved: Some(Unresolved {
law_name: law_name.to_string(),
priority,
point: from,
region: cell_box(config, step, 0, 0, 0),
reason: UnresolvedReason::ReachabilityUndecided {
undecided_cells: undecided,
},
}),
};
};
for (cell, p) in [(a, from), (b, to)] {
if proven_outside[cell.3] {
return Verdict {
violation: Some(Violation {
law_name: law_name.to_string(),
priority,
residual: sdf_eval(node, p),
point: p,
region: cell_box(config, step, cell.0, cell.1, cell.2),
evidence: Evidence::Proved,
}),
unresolved: None,
};
}
}
if flood_reaches(n, a, (b.0, b.1, b.2), &|i| proven_inside[i]) {
return Verdict {
violation: None,
unresolved: None,
};
}
if !flood_reaches(n, a, (b.0, b.1, b.2), &|i| !proven_outside[i]) {
return Verdict {
violation: Some(Violation {
law_name: law_name.to_string(),
priority,
residual: -(from - to).length(),
point: to,
region: cell_box(config, step, b.0, b.1, b.2),
evidence: Evidence::Proved,
}),
unresolved: None,
};
}
Verdict {
violation: None,
unresolved: Some(Unresolved {
law_name: law_name.to_string(),
priority,
point: to,
region: cell_box(config, step, b.0, b.1, b.2),
reason: UnresolvedReason::ReachabilityUndecided {
undecided_cells: undecided,
},
}),
}
}
fn grid_indices(n: usize) -> impl Iterator<Item = (usize, usize, usize, usize)> {
(0..n).flat_map(move |iz| {
(0..n).flat_map(move |iy| (0..n).map(move |ix| (ix, iy, iz, ix + iy * n + iz * n * n)))
})
}
fn check_volume_conservation(
before: &SdfNode,
after: &SdfNode,
relative_tolerance: f32,
law_name: &str,
priority: Priority,
config: &CheckConfig,
) -> Verdict {
let n = config.resolution;
let step = grid_step(config);
let sb = classify_cells(before, config, step);
let sa = classify_cells(after, config, step);
let (mut delta_lo, mut delta_hi) = (0i64, 0i64);
let (mut vol_lo, mut vol_hi) = (0usize, 0usize);
let mut witness: Option<(Vec3, Vec3Interval)> = None;
let frac_bounds = |s: CellSign| match s {
CellSign::Inside => (1i64, 1i64),
CellSign::Outside => (0, 0),
CellSign::Undecided => (0, 1),
};
for (ix, iy, iz, idx) in grid_indices(n) {
let (b_lo, b_hi) = frac_bounds(sb[idx]);
let (a_lo, a_hi) = frac_bounds(sa[idx]);
match sb[idx] {
CellSign::Inside => {
vol_lo += 1;
vol_hi += 1;
}
CellSign::Undecided => vol_hi += 1,
CellSign::Outside => {}
}
delta_lo += a_lo - b_hi;
delta_hi += a_hi - b_lo;
if witness.is_none()
&& matches!(
(sb[idx], sa[idx]),
(CellSign::Outside, CellSign::Inside) | (CellSign::Inside, CellSign::Outside)
)
{
let r = cell_box(config, step, ix, iy, iz);
witness = Some((box_center(r), r));
}
}
if vol_hi == 0 {
return Verdict {
violation: None,
unresolved: None, };
}
let (abs_lo, abs_hi) = if delta_lo <= 0 && delta_hi >= 0 {
(0i64, delta_lo.abs().max(delta_hi))
} else {
let (x, y) = (delta_lo.abs(), delta_hi.abs());
(x.min(y), x.max(y))
};
#[allow(clippy::cast_precision_loss)]
let (diff_lo, diff_hi) = (abs_lo as f32, abs_hi as f32);
#[allow(clippy::cast_precision_loss)]
let (v_lo, v_hi) = (vol_lo as f32, vol_hi as f32);
let rel_lo = diff_lo / v_hi;
let rel_hi = if v_lo > 0.0 {
diff_hi / v_lo
} else {
f32::INFINITY
};
let (point, region) = witness.unwrap_or_else(|| {
let r = Vec3Interval::from_bounds(config.aabb_min, config.aabb_max);
(box_center(r), r)
});
if rel_lo > relative_tolerance {
return Verdict {
violation: Some(Violation {
law_name: law_name.to_string(),
priority,
residual: relative_tolerance - rel_lo, point,
region,
evidence: Evidence::Modelled {
model: MODEL_VOLUME,
},
}),
unresolved: None,
};
}
if rel_hi <= relative_tolerance {
return Verdict {
violation: None,
unresolved: None,
};
}
Verdict {
violation: None,
unresolved: Some(Unresolved {
law_name: law_name.to_string(),
priority,
point,
region,
reason: UnresolvedReason::VolumeUnbracketed {
lo: rel_lo,
hi: rel_hi,
},
}),
}
}
#[derive(Debug, Clone, Default)]
pub struct LawSet {
laws: Vec<Law>,
}
impl LawSet {
#[must_use]
pub fn new() -> Self {
Self::default()
}
#[allow(clippy::missing_const_for_fn)]
pub fn hard(
mut self,
name: impl Into<String>,
constraint: Constraint,
) -> Result<Self, NotProvable> {
self.laws.push(Law::hard(name, constraint)?);
Ok(self)
}
#[must_use]
pub fn soft(mut self, name: impl Into<String>, weight: f32, constraint: Constraint) -> Self {
self.laws.push(Law::soft(name, weight, constraint));
self
}
fn push_provable(mut self, name: impl Into<String>, constraint: Constraint) -> Self {
self.laws.push(Law::hard_unchecked(name, constraint));
self
}
#[must_use]
pub fn laws(&self) -> &[Law] {
&self.laws
}
#[must_use]
pub fn stress(
self,
name: impl Into<String>,
weight: f32,
node: SdfNode,
load_points: Vec<(Vec3, f32)>,
min_thickness_factor: f32,
) -> Self {
self.soft(
name,
weight,
Constraint::Stress {
node,
load_points,
min_thickness_factor,
},
)
}
#[must_use]
pub fn thermal(
self,
name: impl Into<String>,
weight: f32,
node: SdfNode,
heat_sources: Vec<Vec3>,
search_radius: f32,
min_surface_ratio: f32,
) -> Self {
self.soft(
name,
weight,
Constraint::Thermal {
node,
heat_sources,
search_radius,
min_surface_ratio,
},
)
}
#[must_use]
pub fn contact(
self,
name: impl Into<String>,
a: SdfNode,
b: SdfNode,
min_distance: f32,
max_distance: f32,
) -> Self {
self.push_provable(
name,
Constraint::Contact {
a,
b,
min_distance,
max_distance,
},
)
}
#[must_use]
pub fn continuity(
self,
name: impl Into<String>,
weight: f32,
node: SdfNode,
seed_point: Vec3,
) -> Self {
self.soft(name, weight, Constraint::Continuity { node, seed_point })
}
#[must_use]
pub fn gradient_bound(
self,
name: impl Into<String>,
node: SdfNode,
max_gradient: f32,
probe: f32,
) -> Self {
self.push_provable(
name,
Constraint::GradientBound {
node,
max_gradient,
probe,
},
)
}
#[must_use]
pub fn reachable(self, name: impl Into<String>, node: SdfNode, from: Vec3, to: Vec3) -> Self {
self.push_provable(name, Constraint::Reachable { node, from, to })
}
#[must_use]
pub fn volume_conservation(
self,
name: impl Into<String>,
weight: f32,
before: SdfNode,
after: SdfNode,
relative_tolerance: f32,
) -> Self {
self.soft(
name,
weight,
Constraint::VolumeConservation {
before,
after,
relative_tolerance,
},
)
}
#[cfg(feature = "physics")]
#[must_use]
#[allow(clippy::too_many_arguments)]
pub fn thermal_field(
self,
name: impl Into<String>,
weight: f32,
node: SdfNode,
sources: Vec<(Vec3, f32)>,
material: alice_physics::transient_thermal::ThermalMaterial,
metres_per_unit: f32,
ambient_c: f32,
max_temperature_c: f32,
duration_s: f32,
) -> Self {
self.soft(
name,
weight,
Constraint::ThermalField {
node,
sources,
material,
metres_per_unit,
ambient_c,
max_temperature_c,
duration_s,
},
)
}
#[must_use]
pub fn check(&self, config: &CheckConfig) -> LawReport {
check_laws(&self.laws, config)
}
#[must_use]
pub fn detect_contradictions(&self) -> Vec<Contradiction> {
let mut contradictions = Vec::new();
for (i, law_i) in self.laws.iter().enumerate() {
for law_j in &self.laws[i + 1..] {
if let Some(reason) = check_contradiction(&law_i.constraint, &law_j.constraint) {
contradictions.push(Contradiction {
law_a: law_i.name.clone(),
law_b: law_j.name.clone(),
reason,
});
}
}
}
contradictions
}
}
#[derive(Debug, Clone)]
pub struct Contradiction {
pub law_a: String,
pub law_b: String,
pub reason: String,
}
fn check_contradiction(a: &Constraint, b: &Constraint) -> Option<String> {
match (a, b) {
(Constraint::NonOverlap { a: na, b: nb }, Constraint::Containment { inner, outer })
| (Constraint::Containment { inner, outer }, Constraint::NonOverlap { a: na, b: nb }) => {
let dbg_a = format!("{na:?}");
let dbg_b = format!("{nb:?}");
let dbg_inner = format!("{inner:?}");
let dbg_outer = format!("{outer:?}");
if (dbg_a == dbg_inner && dbg_b == dbg_outer)
|| (dbg_a == dbg_outer && dbg_b == dbg_inner)
{
Some(
"NonOverlap と Containment が同一ノードペアに適用: 内包されるならば必ず重なる"
.to_string(),
)
} else {
None
}
}
(
Constraint::NonOverlap { a: na, b: nb },
Constraint::Contact {
a: ca,
b: cb,
min_distance,
..
},
)
| (
Constraint::Contact {
a: ca,
b: cb,
min_distance,
..
},
Constraint::NonOverlap { a: na, b: nb },
) => {
let dbg_a = format!("{na:?}");
let dbg_b = format!("{nb:?}");
let dbg_constraint_a = format!("{ca:?}");
let dbg_constraint_b = format!("{cb:?}");
let same_pair = (dbg_a == dbg_constraint_a && dbg_b == dbg_constraint_b)
|| (dbg_a == dbg_constraint_b && dbg_b == dbg_constraint_a);
if same_pair && *min_distance <= 0.0 {
Some(
"NonOverlap と Contact(min<=0) が同一ノードペアに適用: 接触を許容と禁止が同時"
.to_string(),
)
} else {
None
}
}
_ => None,
}
}
#[must_use]
pub fn top_violations(report: &LawReport, n: usize) -> Vec<&Violation> {
report.violations.iter().take(n).collect()
}
#[must_use]
pub fn hard_violations(report: &LawReport) -> Vec<&Violation> {
report
.violations
.iter()
.filter(|v| v.priority == Priority::Hard)
.collect()
}
#[must_use]
pub fn soft_violations(report: &LawReport) -> Vec<&Violation> {
report
.violations
.iter()
.filter(|v| matches!(v.priority, Priority::Soft(_)))
.collect()
}
#[must_use]
pub fn format_report(report: &LawReport) -> String {
use std::fmt::Write;
let mut out = String::new();
let _ = writeln!(
out,
"Law Check: {}/{} passed",
report.passed, report.total_laws
);
if report.all_passed() {
let _ = writeln!(out, " All laws satisfied (proven at every sample point).");
return out;
}
for u in &report.unresolved {
let _ = writeln!(
out,
" [UNDECIDED] {}: {:?} at=({:.2},{:.2},{:.2}), region=[{:.2}..{:.2}]x[{:.2}..{:.2}]x[{:.2}..{:.2}]",
u.law_name,
u.reason,
u.point.x, u.point.y, u.point.z,
u.region.x.lo, u.region.x.hi,
u.region.y.lo, u.region.y.hi,
u.region.z.lo, u.region.z.hi,
);
}
for v in &report.violations {
let severity = match v.priority {
Priority::Hard => "ERROR",
Priority::Soft(_) => "WARN ",
};
let basis = match v.evidence {
Evidence::Proved => "proved",
Evidence::Modelled { .. } => "modelled",
_ => "witnessed",
};
let _ = writeln!(
out,
" [{severity}] {}: residual={:.4}, basis={basis}, at=({:.2},{:.2},{:.2}), region=[{:.2}..{:.2}]x[{:.2}..{:.2}]x[{:.2}..{:.2}]",
v.law_name,
v.residual,
v.point.x, v.point.y, v.point.z,
v.region.x.lo, v.region.x.hi,
v.region.y.lo, v.region.y.hi,
v.region.z.lo, v.region.z.hi,
);
if let Some(model) = v.evidence.model() {
let _ = writeln!(out, " model: {model}");
}
}
out
}
#[cfg(test)]
mod core_probe_tests {
use super::*;
#[test]
fn an_interval_touching_zero_from_below_is_not_uniformly_negative() {
assert_eq!(interval_sign(Interval { lo: -1.0, hi: 0.0 }), None);
assert_eq!(interval_sign(Interval { lo: -1.0, hi: -0.5 }), Some(false));
assert_eq!(interval_sign(Interval { lo: 0.0, hi: 1.0 }), Some(true));
assert_eq!(interval_sign(Interval::new(0.0, 1.0)), None);
}
#[test]
fn probe_ball_finds_a_surface_inside_the_ball() {
let sphere = SdfNode::sphere(1.0);
let probe = probe_ball(&sphere, Vec3::new(2.0, 0.0, 0.0), 1.5);
match probe {
BallProbe::Crossing { distance } => {
assert!(
(distance - 1.0).abs() < 1e-1,
"表面までの距離は 1.0 のはずが {distance}"
);
}
other => panic!("球内の表面を拾えていない: {other:?}"),
}
}
#[test]
fn probe_ball_returns_the_nearer_crossing_not_the_far_one() {
let left = SdfNode::sphere(1.0).translate(-3.0, 0.0, 0.0);
let right = SdfNode::sphere(1.0).translate(3.0, 0.0, 0.0);
let pair = SdfNode::union(left, right);
let probe = probe_ball(&pair, Vec3::ZERO, 5.0);
match probe {
BallProbe::Crossing { distance } => {
assert!(
(2.0..3.0).contains(&distance),
"近い側 (厳密 2.0、octree 精度で 2.1 前後) を選ぶべきが {distance} \
— 4.0 付近なら遠い交点を採用している"
);
}
other => panic!("交点を拾えていない: {other:?}"),
}
}
#[test]
fn gap_exceeds_is_true_and_subdivides_for_a_gap_just_above_m() {
let a = SdfNode::sphere(1.0).translate(-2.0, 0.0, 0.0);
let b = SdfNode::sphere(1.0).translate(2.0, 0.0, 0.0);
let config = CheckConfig {
aabb_min: Vec3::splat(-3.0),
aabb_max: Vec3::splat(3.0),
resolution: 3,
};
assert!(
gap_exceeds(&a, &b, 1.5, &config),
"gap 2.0 の 2 球で m = 1.5 が満たされない"
);
}
#[test]
fn gap_exceeds_does_not_clear_a_box_whose_interval_touches_zero() {
let boxy = SdfNode::box3d(2.0, 2.0, 2.0);
let everywhere = SdfNode::sphere(50.0);
let config = CheckConfig {
aabb_min: Vec3::splat(1.5),
aabb_max: Vec3::splat(2.0),
resolution: 1,
};
assert!(
!gap_exceeds(&boxy, &everywhere, 1.0, &config),
"a 側の区間下界 0.0 (表面に接する) を「離れている」と刈っている"
);
assert!(
!gap_exceeds(&everywhere, &boxy, 1.0, &config),
"b 側の区間下界 0.0 (表面に接する) を「離れている」と刈っている"
);
}
#[test]
fn contact_upper_bound_radius_includes_the_cell_diagonal() {
let a = SdfNode::sphere(1.0);
let b = SdfNode::sphere(1.0).translate(4.0, 0.0, 0.0);
let config = CheckConfig {
aabb_min: Vec3::splat(-2.0),
aabb_max: Vec3::splat(2.0),
resolution: 2,
};
let found = contact_upper_bound(&a, &b, 0.1, 0.2, &config);
let Some((ub, point, _)) = found else {
panic!("cell 対角を含む半径なら上界が取れるはずが None");
};
assert!(ub >= 2.0, "上界 {ub} が真の gap 2.0 を下回っている");
assert!(
sdf_eval(&a, point) > 0.0 && sdf_eval(&b, point) > 0.0,
"上界の標本点 {point:?} が形状の内部にある"
);
}
#[test]
fn contact_upper_bound_only_samples_points_outside_both() {
let a = SdfNode::sphere(1.0);
let b = SdfNode::sphere(1.0).translate(2.6, 0.0, 0.0);
let config = CheckConfig {
aabb_min: Vec3::splat(-1.0),
aabb_max: Vec3::splat(1.0),
resolution: 2,
};
assert!(
contact_upper_bound(&a, &b, 0.1, 0.2, &config).is_none(),
"標本点が全て a の内部なのに上界を返した"
);
}
#[allow(clippy::float_cmp)]
#[test]
fn a_point_on_either_boundary_is_not_a_containment_witness() {
const EPS: f32 = 1.0e-3;
let config = CheckConfig {
aabb_min: Vec3::new(0.5, -0.5, -0.5),
aabb_max: Vec3::new(1.5, 0.5, 0.5),
resolution: 1,
};
let boundary = Vec3::new(1.0, 0.0, 0.0);
let inner = SdfNode::box3d(2.0, 2.0, 2.0);
let outer = SdfNode::box3d(2.0f32.mul_add(-EPS, 2.0), 2.0, 2.0);
assert_eq!(
sdf_eval(&inner, boundary),
0.0,
"inner の表面が厳密 0 でない"
);
assert!(
sdf_eval(&outer, boundary) > 0.0,
"outer の外部になっていない"
);
let report = check_laws(
&[
Law::hard("surface_of_inner", Constraint::Containment { inner, outer })
.expect("provable constraint"),
],
&config,
);
assert!(
!report.has_hard_violations(),
"inner の表面上の点 (fi = 0) をはみ出しの witness にした\n{}",
format_report(&report)
);
let inner2 = SdfNode::box3d(2.0f32.mul_add(EPS, 2.0), 2.0, 2.0);
let outer2 = SdfNode::box3d(2.0, 2.0, 2.0);
assert!(
sdf_eval(&inner2, boundary) < 0.0,
"inner の内部になっていない"
);
assert_eq!(
sdf_eval(&outer2, boundary),
0.0,
"outer の表面が厳密 0 でない"
);
let report2 = check_laws(
&[Law::hard(
"surface_of_outer",
Constraint::Containment {
inner: inner2,
outer: outer2,
},
)
.expect("provable constraint")],
&config,
);
assert!(
!report2.has_hard_violations(),
"outer の表面上の点 (fo = 0) をはみ出しの witness にした\n{}",
format_report(&report2)
);
}
#[allow(clippy::float_cmp)]
#[test]
fn contact_upper_bound_equal_to_min_distance_is_not_too_close() {
let a = SdfNode::box3d(2.0, 2.0, 2.0);
let b = SdfNode::box3d(2.0, 2.0, 2.0).translate(3.0, 0.0, 0.0);
let config = CheckConfig {
aabb_min: Vec3::new(1.0, -0.5, -0.5),
aabb_max: Vec3::new(2.0, 0.5, 0.5),
resolution: 1,
};
let mid = Vec3::new(1.5, 0.0, 0.0);
assert_eq!(sdf_eval(&a, mid), 0.5, "a までの距離が 0.5 でない");
assert_eq!(sdf_eval(&b, mid), 0.5, "b までの距離が 0.5 でない");
let report = check_laws(
&[Law::hard(
"gap_at_the_lower_bound",
Constraint::Contact {
a,
b,
min_distance: 1.0,
max_distance: 2.0,
},
)
.expect("provable constraint")],
&config,
);
assert!(
!report.has_hard_violations(),
"gap が下限ちょうど (1.0) の scene を「近すぎ」と報告した\n{}",
format_report(&report)
);
}
#[test]
fn gap_exceeds_is_false_for_touching_shapes() {
let a = SdfNode::sphere(1.0).translate(-1.0, 0.0, 0.0);
let b = SdfNode::sphere(1.0).translate(1.0, 0.0, 0.0);
let config = CheckConfig::default();
assert!(
!gap_exceeds(&a, &b, 0.5, &config),
"原点で接している 2 球が gap 0.5 を満たしてしまっている"
);
}
#[allow(clippy::float_cmp)]
#[test]
fn a_point_on_one_surface_is_not_an_overlap_witness() {
const EPS: f32 = 1.0e-3;
let config = CheckConfig {
aabb_min: Vec3::new(0.5, -0.5, -0.5),
aabb_max: Vec3::new(1.5, 0.5, 0.5),
resolution: 1,
};
let touch = Vec3::new(1.0, 0.0, 0.0);
let a = SdfNode::box3d(2.0, 2.0, 2.0);
let b = SdfNode::box3d(2.0, 2.0, 2.0).translate(2.0 - EPS, 0.0, 0.0);
assert_eq!(sdf_eval(&a, touch), 0.0, "a の表面の場が厳密 0 でない");
assert!(sdf_eval(&b, touch) < 0.0, "b の内部になっていない");
let report = check_laws(
&[Law::hard("surface_of_a", Constraint::NonOverlap { a, b })
.expect("provable constraint")],
&config,
);
assert!(
!report.has_hard_violations(),
"a の表面上の点 (fa = 0) を侵入の witness にした\n{}",
format_report(&report)
);
let a2 = SdfNode::box3d(2.0f32.mul_add(EPS, 2.0), 2.0, 2.0);
let b2 = SdfNode::box3d(2.0, 2.0, 2.0).translate(2.0, 0.0, 0.0);
assert_eq!(sdf_eval(&b2, touch), 0.0, "b の表面の場が厳密 0 でない");
assert!(sdf_eval(&a2, touch) < 0.0, "a の内部になっていない");
let report2 = check_laws(
&[
Law::hard("surface_of_b", Constraint::NonOverlap { a: a2, b: b2 })
.expect("provable constraint"),
],
&config,
);
assert!(
!report2.has_hard_violations(),
"b の表面上の点 (fb = 0) を侵入の witness にした\n{}",
format_report(&report2)
);
}
#[test]
fn the_lipschitz_enclosure_uses_the_half_diagonal_and_the_product_with_l() {
let bx = Vec3Interval {
x: Interval { lo: 0.5, hi: 2.5 },
y: Interval { lo: 1.0, hi: 3.0 },
z: Interval { lo: 0.25, hi: 1.25 },
};
let e = lipschitz_enclosure(bx, 10.0, 2.0)
.expect("外部 (f(c) > 0) で上界も有限なので包囲が出る");
assert!(
(e.lo - 7.0).abs() < 1e-3 && (e.hi - 13.0).abs() < 1e-3,
"半対角 1.5 × L 2.0 = 余裕 3.0 で [7, 13] のはずが [{}, {}]",
e.lo,
e.hi
);
}
#[test]
fn the_lipschitz_enclosure_refuses_interior_boxes_and_useless_bounds() {
let bx = Vec3Interval {
x: Interval { lo: 0.5, hi: 2.5 },
y: Interval { lo: 1.0, hi: 3.0 },
z: Interval { lo: 0.25, hi: 1.25 },
};
assert!(
lipschitz_enclosure(bx, -0.1, 2.0).is_none(),
"f(c) < 0 で主張した"
);
assert!(
lipschitz_enclosure(bx, 10.0, 0.0).is_none(),
"L = 0 で主張した"
);
assert!(
lipschitz_enclosure(bx, 10.0, -1.0).is_none(),
"L < 0 で主張した"
);
assert!(
lipschitz_enclosure(bx, 10.0, f32::INFINITY).is_none(),
"L が非有限で主張した"
);
assert!(
lipschitz_enclosure(bx, f32::NAN, 2.0).is_none(),
"f(c) が NaN で主張した"
);
assert!(
lipschitz_enclosure(bx, 0.0, 2.0).is_some(),
"f(c) = 0 は契約内 (`f < 0` で弾く) なのに主張しなかった"
);
}
#[test]
fn a_single_point_intersection_is_the_tightest_enclosure_not_a_contradiction() {
let bx = Vec3Interval {
x: Interval { lo: 0.5, hi: 2.5 },
y: Interval { lo: 1.0, hi: 3.0 },
z: Interval { lo: 0.25, hi: 1.25 },
};
let iv = Interval { lo: 10.0, hi: 10.0 };
match refine_with_lipschitz(bx, iv, 10.0, 2.0) {
Refined::Enclosure(r) => assert!(
(r.lo - 10.0).abs() < 1e-3 && (r.hi - 10.0).abs() < 1e-3,
"1 点交差が [{}, {}] になった",
r.lo,
r.hi
),
Refined::Contradiction => panic!("1 点で接する交差を矛盾として扱った"),
}
}
#[test]
fn an_empty_intersection_is_reported_as_a_contradiction() {
let bx = Vec3Interval {
x: Interval { lo: 0.5, hi: 2.5 },
y: Interval { lo: 1.0, hi: 3.0 },
z: Interval { lo: 0.25, hi: 1.25 },
};
let iv = Interval {
lo: 100.0,
hi: 200.0,
};
assert!(
matches!(
refine_with_lipschitz(bx, iv, 10.0, 2.0),
Refined::Contradiction
),
"重ならない 2 包囲を矛盾として扱わなかった"
);
}
#[test]
fn without_a_lipschitz_bound_the_interval_passes_through_unchanged() {
let bx = Vec3Interval {
x: Interval { lo: 0.5, hi: 2.5 },
y: Interval { lo: 1.0, hi: 3.0 },
z: Interval { lo: 0.25, hi: 1.25 },
};
let iv = Interval { lo: -5.0, hi: 5.0 };
match refine_with_lipschitz(bx, iv, -0.1, 2.0) {
Refined::Enclosure(r) => {
assert!(
(r.lo - iv.lo).abs() < 1e-6 && (r.hi - iv.hi).abs() < 1e-6,
"契約外なのに区間が変わった [{}, {}]",
r.lo,
r.hi
);
}
Refined::Contradiction => panic!("締める材料が無い時に矛盾を報告した"),
}
}
}
#[cfg(test)]
mod evidence_gate_tests {
use super::*;
fn every_constraint() -> Vec<Constraint> {
let n = || SdfNode::sphere(1.0);
vec![
Constraint::NonOverlap { a: n(), b: n() },
Constraint::Containment {
inner: n(),
outer: n(),
},
Constraint::MinThickness {
node: n(),
min_thickness: 0.5,
},
Constraint::Stress {
node: n(),
load_points: vec![(Vec3::ZERO, 1.0)],
min_thickness_factor: 0.2,
},
Constraint::Thermal {
node: n(),
heat_sources: vec![Vec3::ZERO],
search_radius: 1.0,
min_surface_ratio: 0.8,
},
Constraint::Contact {
a: n(),
b: n(),
min_distance: 0.1,
max_distance: 1.0,
},
Constraint::Continuity {
node: n(),
seed_point: Vec3::ZERO,
},
Constraint::GradientBound {
node: n(),
max_gradient: 1.0,
probe: 0.01,
},
Constraint::Reachable {
node: n(),
from: Vec3::ZERO,
to: Vec3::ZERO,
},
Constraint::VolumeConservation {
before: n(),
after: n(),
relative_tolerance: 0.05,
},
]
}
#[test]
fn all_variants_are_classified() {
assert_eq!(
every_constraint().len(),
10,
"variant を足したら every_constraint にも足す"
);
let modelled: Vec<&str> = every_constraint()
.iter()
.filter(|c| !c.evidence_class().is_proof())
.map(Constraint::name)
.collect();
assert_eq!(
modelled,
["Stress", "Thermal", "Continuity", "VolumeConservation"],
"モデル推定の顔ぶれが変わった — 意図した変更なら CHANGELOG に書く"
);
}
#[test]
fn hard_is_gated_by_evidence_class() {
for c in every_constraint() {
let (name, class) = (c.name(), c.evidence_class());
let got = Law::hard("x", c);
assert_eq!(
got.is_ok(),
class.is_proof(),
"{name} の gate が evidence_class と食い違う"
);
if let Err(e) = got {
assert_eq!(e.constraint, name);
assert_eq!(Some(e.model), class.model(), "{name} の model 説明が違う");
assert!(
e.to_string().contains("Law::soft"),
"{name}: 対処法が message に無い"
);
}
}
}
#[test]
fn modelled_constraints_cannot_produce_hard_violations() {
let split = SdfNode::sphere(0.8).union(SdfNode::sphere(0.8).translate(4.0, 0.0, 0.0));
let config = CheckConfig {
aabb_min: Vec3::new(-2.0, -2.0, -2.0),
aabb_max: Vec3::new(6.0, 2.0, 2.0),
resolution: 16,
};
let report = LawSet::new()
.continuity("split", 1.0, split.clone(), Vec3::ZERO)
.volume_conservation("shrunk", 1.0, split, SdfNode::sphere(0.8), 0.05)
.check(&config);
assert!(
!report.violations.is_empty(),
"検査対象が違反を出していない (test が無意味になっている)"
);
assert!(
!report.has_hard_violations(),
"モデル推定が Hard 違反として出た\n{}",
format_report(&report)
);
assert!(
report.proven_violations().is_empty(),
"モデル推定が proven_violations に混ざった"
);
for v in &report.violations {
assert!(
v.evidence.model().is_some(),
"{} の違反が Modelled でない",
v.law_name
);
}
}
#[test]
fn reported_evidence_never_exceeds_the_declared_class() {
let config = CheckConfig {
aabb_min: Vec3::splat(-2.0),
aabb_max: Vec3::splat(2.0),
resolution: 8,
};
let overlap = SdfNode::sphere(1.0);
let cases: Vec<Law> = vec![
Law::hard(
"overlap",
Constraint::NonOverlap {
a: overlap.clone(),
b: overlap.clone(),
},
)
.expect("provable constraint"),
Law::soft(
"heat",
1.0,
Constraint::Thermal {
node: overlap.clone(),
heat_sources: vec![Vec3::ZERO],
search_radius: 1.5,
min_surface_ratio: 0.99,
},
),
Law::soft(
"load",
1.0,
Constraint::Stress {
node: overlap,
load_points: vec![(Vec3::ZERO, 5.0)],
min_thickness_factor: 0.5,
},
),
];
let classes: Vec<(String, Evidence)> = cases
.iter()
.map(|l| (l.name.clone(), l.constraint.evidence_class()))
.collect();
let report = check_laws(&cases, &config);
assert!(
!report.violations.is_empty(),
"どの法則も違反を出さなかった"
);
for v in &report.violations {
let (_, class) = classes
.iter()
.find(|(n, _)| *n == v.law_name)
.expect("報告された法則名が入力に無い");
assert!(
!v.evidence.is_proof() || class.is_proof(),
"{}: 実際の根拠 {:?} が申告した class {:?} より強い — \
evidence_class の表が実装から遅れている",
v.law_name,
v.evidence,
class
);
}
}
#[test]
fn surface_band_does_not_depend_on_the_field_scale() {
let node = SdfNode::gyroid(0.6, 1.2);
let config = CheckConfig {
aabb_min: Vec3::splat(-2.0),
aabb_max: Vec3::splat(2.0),
resolution: 24,
};
let (src, radius) = (Vec3::ZERO, 1.5_f32);
let step = grid_step(&config);
let signs = classify_cells(&node, &config, step);
let (mut straddling, mut missed_by_field, mut inside) = (0usize, 0usize, 0usize);
for (ix, iy, iz, idx) in grid_indices(config.resolution) {
let center = cell_center(&config, step, ix, iy, iz);
if center.distance(src) > radius {
continue;
}
if signs[idx] == CellSign::Inside {
inside += 1;
}
let region = cell_box(&config, step, ix, iy, iz);
if !cell_straddles_surface(&node, region) {
continue;
}
straddling += 1;
assert_eq!(
signs[idx],
CellSign::Undecided,
"表面を含むセルが内外確定になった (帯が上界でない)"
);
if sdf_eval(&node, center).abs() >= step.x {
missed_by_field += 1;
}
}
assert!(
straddling > 0 && inside > 0,
"検査対象が空 (straddling={straddling} inside={inside})"
);
assert!(
missed_by_field > 0,
"場の帯と区間の帯が一致してしまい基準の差を測れない — \
場が過大申告する node に変える"
);
let (_, hi, _) = surface_ratio_bracket(&node, src, radius, &config, &signs)
.expect("内部確定セルがあるので bracket が取れる");
#[allow(clippy::cast_precision_loss)]
let needed = straddling as f32 / inside as f32;
assert!(
hi >= needed,
"表面比の上界 {hi:.4} が、表面を含むと確定したセル分 {needed:.4} に \
届いていない = 帯が場の値で痩せている \
(straddling={straddling} inside={inside} missed_by_field={missed_by_field})"
);
}
#[cfg(feature = "physics")]
fn constant_material(
k: f32,
cp: f32,
rho: f32,
) -> alice_physics::transient_thermal::ThermalMaterial {
use alice_physics::transient_thermal::{TemperatureDependence as Dep, ThermalMaterial};
ThermalMaterial {
name: "oracle_constant",
conductivity: Dep::Constant(k),
specific_heat: Dep::Constant(cp),
density: Dep::Constant(rho),
reference_temperature: 293.15,
}
}
#[cfg(feature = "physics")]
fn filled_cube(n: usize, half: f32) -> (CheckConfig, Vec<bool>) {
let config = CheckConfig {
aabb_min: Vec3::splat(-half),
aabb_max: Vec3::splat(half),
resolution: n,
};
(config, vec![true; n * n * n])
}
#[cfg(feature = "physics")]
#[test]
fn adiabatic_peak_is_bounded_by_the_energy_balance() {
let (k, cp, rho) = (200.0_f32, 900.0_f32, 2700.0_f32); let material = constant_material(k, cp, rho);
let n = 9;
let half = 0.5_f32; let (config, solid) = filled_cube(n, half);
let (watts, ambient) = (50.0_f32, 20.0_f32);
let volume = (2.0 * half).powi(3);
#[allow(clippy::cast_precision_loss)]
let cell_volume = volume / (n * n * n) as f32;
let peak_of = |duration: f32| {
let scene = ThermalScene {
solid: &solid,
sources: &[(Vec3::ZERO, watts)],
material: &material,
metres_per_unit: 1.0,
ambient_c: ambient,
duration_s: duration,
};
solve_peak_temperature(&scene, SurfaceCooling::Adiabatic, &config).expect("解けるはず")
};
let mut previous_ratio = f32::INFINITY;
let mut excess = Vec::new();
for duration in [200.0_f32, 1000.0, 5000.0] {
let peak = peak_of(duration);
let mean_rise = watts * duration / (rho * cp * volume);
let single_cell_rise = watts * duration / (rho * cp * cell_volume);
assert!(
peak >= ambient + mean_rise * 0.999,
"t={duration}: 最高温度 {peak:.4} が平均 {:.4} 未満 = 熱が消えている",
ambient + mean_rise
);
assert!(
peak <= ambient + single_cell_rise,
"t={duration}: 最高温度 {peak:.4} が 1 セル集中の上限 {:.4} 超 = 熱が湧いている",
ambient + single_cell_rise
);
let ratio = (peak - ambient) / mean_rise;
assert!(
ratio < previous_ratio,
"t={duration}: 最高/平均 の比 {ratio:.3} が前回 {previous_ratio:.3} から \
下がらない = 拡散が進んでいない"
);
previous_ratio = ratio;
excess.push(peak - ambient - mean_rise);
}
let (d1, d2) = ((excess[1] - excess[0]).abs(), (excess[2] - excess[1]).abs());
assert!(
d2 * 5.0 < d1,
"超過温度が収束しない (差 {d1:.6} → {d2:.6}) = 定常状態に入っていない"
);
#[allow(clippy::cast_precision_loss)]
let dx = 2.0 * half / n as f32;
let continuum = watts / (6.0 * k * dx);
let rel = excess[2] / continuum;
assert!(
(0.6..=1.4).contains(&rel),
"定常超過 {:.4} が連続体の見積もり {continuum:.4} の {rel:.3} 倍 = \
熱伝導率か単位の扱いがずれている",
excess[2]
);
}
#[cfg(feature = "physics")]
#[test]
fn no_source_keeps_the_field_at_ambient() {
let material = constant_material(200.0, 900.0, 2700.0);
let (config, solid) = filled_cube(9, 0.5);
for cooling in [SurfaceCooling::Adiabatic, SurfaceCooling::Isothermal] {
let scene = ThermalScene {
solid: &solid,
sources: &[],
material: &material,
metres_per_unit: 1.0,
ambient_c: 20.0,
duration_s: 100.0,
};
let peak = solve_peak_temperature(&scene, cooling, &config).expect("解けるはず");
assert!(
(peak - 20.0).abs() < 1.0e-3,
"熱源なしで温度が {peak:.6} に動いた"
);
}
}
#[cfg(feature = "physics")]
#[test]
fn isothermal_never_exceeds_adiabatic() {
let material = constant_material(200.0, 900.0, 2700.0);
let n = 9;
let config = CheckConfig {
aabb_min: Vec3::splat(-0.5),
aabb_max: Vec3::splat(0.5),
resolution: n,
};
let signs = classify_cells(&SdfNode::sphere(0.35), &config, grid_step(&config));
let solid: Vec<bool> = signs.iter().map(|s| *s != CellSign::Outside).collect();
assert!(solid.iter().any(|b| *b), "材料セルが無い");
assert!(
solid.iter().any(|b| !*b),
"外側セルが無い (境界条件が効かない)"
);
let solve = |c| {
let scene = ThermalScene {
solid: &solid,
sources: &[(Vec3::ZERO, 50.0)],
material: &material,
metres_per_unit: 1.0,
ambient_c: 20.0,
duration_s: 200.0,
};
solve_peak_temperature(&scene, c, &config).expect("解けるはず")
};
let (hot, cold) = (
solve(SurfaceCooling::Adiabatic),
solve(SurfaceCooling::Isothermal),
);
assert!(
cold <= hot,
"等温 (h = ∞) の最高温度 {cold:.4} が断熱 (h = 0) の {hot:.4} を超えた \
= 上下界が逆転している"
);
assert!(
cold < hot,
"境界条件を変えても最高温度が動かない ({cold:.4}) = 表面が効いていない"
);
}
#[cfg(feature = "physics")]
fn sphere_in_box(n: usize) -> (CheckConfig, Vec<bool>) {
let config = CheckConfig {
aabb_min: Vec3::splat(-0.5),
aabb_max: Vec3::splat(0.5),
resolution: n,
};
let signs = classify_cells(&SdfNode::sphere(0.35), &config, grid_step(&config));
let solid: Vec<bool> = signs.iter().map(|s| *s != CellSign::Outside).collect();
assert!(solid.iter().any(|b| *b), "材料セルが無い");
assert!(solid.iter().any(|b| !*b), "外側セルが無い");
(config, solid)
}
#[cfg(feature = "physics")]
#[test]
fn isothermal_saturates_but_adiabatic_keeps_rising() {
let material = constant_material(200.0, 900.0, 2700.0);
let (config, solid) = sphere_in_box(9);
let rise = |cooling, duration| {
let scene = ThermalScene {
solid: &solid,
sources: &[(Vec3::ZERO, 50.0)],
material: &material,
metres_per_unit: 1.0,
ambient_c: 20.0,
duration_s: duration,
};
solve_peak_temperature(&scene, cooling, &config).expect("解けるはず") - 20.0
};
let iso_short = rise(SurfaceCooling::Isothermal, 3000.0);
let iso_long = rise(SurfaceCooling::Isothermal, 30000.0);
let adi_short = rise(SurfaceCooling::Adiabatic, 3000.0);
let adi_long = rise(SurfaceCooling::Adiabatic, 30000.0);
assert!(
iso_long < iso_short * 1.2,
"等温境界で時間を 10 倍にしたら温度上昇が {iso_short:.4} → {iso_long:.4} \
(1.2 倍超) = 定常に入っていない = 表面が周囲温度に固定されていない"
);
assert!(
adi_long > adi_short * 3.0,
"断熱境界で時間を 10 倍にしても温度上昇が {adi_short:.4} → {adi_long:.4} \
(3 倍未満) = 熱が逃げている"
);
}
#[cfg(feature = "physics")]
#[test]
fn material_moves_the_peak_in_opposite_directions_per_boundary() {
let material = constant_material(200.0, 900.0, 2700.0);
let config = CheckConfig {
aabb_min: Vec3::splat(-0.5),
aabb_max: Vec3::splat(0.5),
resolution: 9,
};
let signs = classify_cells(&SdfNode::sphere(0.35), &config, grid_step(&config));
let strict: Vec<bool> = signs.iter().map(|s| *s == CellSign::Inside).collect();
let loose: Vec<bool> = signs.iter().map(|s| *s != CellSign::Outside).collect();
let extra = (0..signs.len()).filter(|&i| loose[i] && !strict[i]).count();
assert!(
strict.iter().any(|b| *b) && extra > 0,
"内部確定セルと未確定セルが両方要る (strict={} extra={extra})",
strict.iter().filter(|b| **b).count()
);
let sources = [(Vec3::ZERO, 50.0)];
let peak = |solid: &[bool], cooling| {
let scene = ThermalScene {
solid,
sources: &sources,
material: &material,
metres_per_unit: 1.0,
ambient_c: 20.0,
duration_s: 2000.0,
};
solve_peak_temperature(&scene, cooling, &config).expect("解けるはず")
};
let adi_small = peak(&strict, SurfaceCooling::Adiabatic);
let adi_big = peak(&loose, SurfaceCooling::Adiabatic);
assert!(
adi_big < adi_small,
"断熱で材料を増やしたのに最高温度が下がらない ({adi_small:.4} → {adi_big:.4}) \
= 上界を最小材料に置く根拠が崩れる"
);
let iso_small = peak(&strict, SurfaceCooling::Isothermal);
let iso_big = peak(&loose, SurfaceCooling::Isothermal);
assert!(
iso_big > iso_small,
"等温で材料を増やしたのに最高温度が上がらない ({iso_small:.4} → {iso_big:.4}) \
= 下界を最小材料に置く根拠が崩れる"
);
}
#[cfg(feature = "physics")]
#[test]
fn bracket_contains_every_admissible_material_set() {
let material = constant_material(200.0, 900.0, 2700.0);
let config = CheckConfig {
aabb_min: Vec3::splat(-0.5),
aabb_max: Vec3::splat(0.5),
resolution: 9,
};
let signs = classify_cells(&SdfNode::sphere(0.35), &config, grid_step(&config));
let strict: Vec<bool> = signs.iter().map(|s| *s == CellSign::Inside).collect();
let undecided: Vec<usize> = (0..signs.len())
.filter(|&i| signs[i] == CellSign::Undecided)
.collect();
assert!(!undecided.is_empty(), "未確定セルが無いと挟む意味が無い");
let sources = [(Vec3::ZERO, 50.0)];
let peak = |solid: &[bool], cooling| {
let scene = ThermalScene {
solid,
sources: &sources,
material: &material,
metres_per_unit: 1.0,
ambient_c: 20.0,
duration_s: 2000.0,
};
solve_peak_temperature(&scene, cooling, &config).expect("解けるはず")
};
let probe = peak(&strict, SurfaceCooling::Adiabatic)
.midpoint(peak(&strict, SurfaceCooling::Isothermal));
let report = check_laws(
&[Law::soft(
"heat",
1.0,
Constraint::ThermalField {
node: SdfNode::sphere(0.35),
sources: vec![(Vec3::ZERO, 50.0)],
material,
metres_per_unit: 1.0,
ambient_c: 20.0,
max_temperature_c: probe,
duration_s: 2000.0,
},
)],
&config,
);
let (lo, hi) = match report.unresolved.first().map(|u| &u.reason) {
Some(UnresolvedReason::TemperatureUnbracketed { lo_c, hi_c }) => (*lo_c, *hi_c),
other => panic!("bracket の中点なら未定になるはず (実際 {other:?} / {report:?})"),
};
assert!(lo < hi, "bracket が潰れている ({lo:.4}, {hi:.4})");
for numerator in [0_usize, 1, 2, 3, 4] {
let mut m = strict.clone();
for (rank, &idx) in undecided.iter().enumerate() {
if rank * 4 < numerator * undecided.len() {
m[idx] = true;
}
}
for cooling in [SurfaceCooling::Adiabatic, SurfaceCooling::Isothermal] {
let t = peak(&m, cooling);
assert!(
(lo..=hi).contains(&t),
"未確定セル {numerator}/4 採用 / {cooling:?} の最高温度 {t:.4} が \
bracket [{lo:.4}, {hi:.4}] の外 = 挟めていない"
);
}
}
}
#[cfg(feature = "physics")]
#[test]
fn thermal_solve_is_bit_reproducible() {
let material = constant_material(200.0, 900.0, 2700.0);
let (config, solid) = filled_cube(9, 0.5);
let run = || {
let scene = ThermalScene {
solid: &solid,
sources: &[(Vec3::new(0.1, -0.2, 0.05), 50.0)],
material: &material,
metres_per_unit: 1.0,
ambient_c: 20.0,
duration_s: 150.0,
};
solve_peak_temperature(&scene, SurfaceCooling::Adiabatic, &config).expect("解けるはず")
};
assert_eq!(run().to_bits(), run().to_bits(), "同じ入力で結果が揺れた");
}
fn l_arms() -> (SdfNode, SdfNode) {
(
SdfNode::box3d_half_extents(0.5, 0.15, 1.0).translate(-0.5, -0.15, 0.0),
SdfNode::box3d_half_extents(0.15, 0.5, 1.0).translate(-0.15, -0.5, 0.0),
)
}
#[test]
fn solid_angle_fraction_is_scale_invariant_at_a_sharp_corner() {
let (arm_x, arm_y) = l_arms();
let node = arm_x.union(arm_y);
let corner = Vec3::new(-0.3, -0.3, 0.0);
let seen: Vec<f32> = [0.26_f32, 0.13, 0.065, 0.032]
.iter()
.map(|&r| material_fraction(&node, corner, r))
.collect();
for f in &seen {
assert!(
(0.72..=0.78).contains(f),
"占有率 {f:.3} が内角 3π/2 の α/2π = 0.75 近傍にない: {seen:?}"
);
}
let spread = seen.iter().copied().fold(f32::MIN, f32::max)
- seen.iter().copied().fold(f32::MAX, f32::min);
assert!(
spread < 0.02,
"尺度を 8 倍変えて占有率が {spread:.3} 動いた = 尺度不変でない: {seen:?}"
);
}
#[test]
fn solid_angle_fraction_decays_with_radius_on_a_fillet() {
let (arm_x, arm_y) = l_arms();
let node = arm_x.smooth_union(arm_y, 0.35);
let mut surf = Vec3::new(-0.45, -0.45, 0.0);
for _ in 0..12 {
let value = sdf_eval(&node, surf);
let grad = alice_sdf::eval::eval_gradient(&node, surf);
let norm_sq = grad.length_squared();
assert!(norm_sq > 1.0e-12, "勾配が消えた = 射影できない");
surf -= grad * (value / norm_sq);
}
assert!(
sdf_eval(&node, surf).abs() < 1.0e-3,
"フィレット面へ射影できていない (f = {})",
sdf_eval(&node, surf)
);
for radius in [0.26_f32, 0.173, 0.13] {
let d_far = material_fraction(&node, surf, radius) - 0.5;
let d_near = material_fraction(&node, surf, radius * 0.5) - 0.5;
assert!(
d_far > 0.01,
"r={radius}: フィレット面なのに隔たり {d_far:.3} が小さすぎて比を測れない"
);
let ratio = d_near / d_far;
assert!(
(0.35..=0.75).contains(&ratio),
"r={radius}: 比 {ratio:.3} が理論値 0.5 の近傍にない \
(鋭い入隅なら 1.0 になるので、0.75 以上だと区別が付かない)"
);
}
}
#[test]
fn lawset_convenience_keeps_modelled_laws_soft() {
let n = || SdfNode::sphere(1.0);
let set = LawSet::new()
.stress("s", 0.5, n(), vec![(Vec3::ZERO, 1.0)], 0.2)
.thermal("t", 0.5, n(), vec![Vec3::ZERO], 1.0, 0.8)
.continuity("c", 0.5, n(), Vec3::ZERO)
.volume_conservation("v", 0.5, n(), n(), 0.05);
assert_eq!(set.laws().len(), 4);
for law in set.laws() {
assert_eq!(
law.priority,
Priority::Soft(0.5),
"{} が Hard で積まれた",
law.name
);
}
let provable = LawSet::new()
.contact("k", n(), n(), 0.1, 1.0)
.gradient_bound("g", n(), 1.0, 0.01)
.reachable("r", n(), Vec3::ZERO, Vec3::ZERO);
for law in provable.laws() {
assert_eq!(law.priority, Priority::Hard, "{} が降格した", law.name);
}
}
}