Merge pull request #1595 from brett19/port/brep-transform-matrix

port: SmBrep::Transform with a 4x4 matrix and SmMatrix::SVDecompose (#1570)
This commit is contained in:
bherbruck
2026-10-30 17:05:08 -04:00
committed by GitHub
13 changed files with 2966 additions and 18 deletions
+102 -2
View File
@@ -35,7 +35,7 @@ const SM_ERR_WARNING: i32 = 1023;
pub struct Runner;
/// Operation names used in `inputs.op`.
const OPS: [(&str, Op); 11] = [
const OPS: [(&str, Op); 12] = [
("copy", Op::Copy),
("resize", Op::Resize),
("access", Op::Access),
@@ -47,8 +47,13 @@ const OPS: [(&str, Op); 11] = [
("solve", Op::Solve),
("eigen", Op::Eigen),
("eigen3", Op::Eigen3),
("svd", Op::Svd),
];
/// Generated `SVDecompose` cases per run (named `svd####`, seeded apart from the
/// `gen####` stream so its cases keep their inputs).
pub const SVD_GENERATED: usize = 600;
#[derive(Debug, Clone, Serialize, Deserialize)]
struct Inputs {
op: String,
@@ -219,6 +224,13 @@ fn rust_run(op: Op, i: &Inputs) -> Result<solidcore::Result<Output>, HarnessErro
out.extend_from_slice(a.as_slice());
result(out, Vec::new(), None)
}
Op::Svd => {
let (w, v) = a.sv_decompose()?;
let mut out = w;
out.extend_from_slice(v.as_slice());
out.extend_from_slice(a.as_slice());
result(out, Vec::new(), None)
}
})
})())
}
@@ -473,6 +485,82 @@ fn gen_case(rng: &mut Rng) -> Inputs {
}
}
/// An `SVDecompose` input: mostly `m >= n` (sometimes `m < n`, `SM_ERR`), with the
/// shapes of [`square`], 3x3 transform blocks (rotation, shear, mirror, scale,
/// near-singular) and huge / tiny scales. Entries the FFI guard refuses become zero.
fn gen_svd(rng: &mut Rng) -> Inputs {
let (rows, cols, mut m) = if rng.chance(0.3) {
(3, 3, transform_block(rng))
} else {
let rows = dim(rng);
let cols = if rng.chance(0.05) {
rows + 1
} else if rng.chance(0.5) {
rows
} else {
1 + rng.below(rows as u64) as usize
};
let m = if rows == cols {
let symmetric = rng.chance(0.2);
square(rng, rows, symmetric)
} else {
(0..rows * cols).map(|_| entry(rng)).collect()
};
(rows, cols, m)
};
if rng.chance(0.1) {
let k = 10f64.powi(rng.range_i64(-60, 60) as i32);
for x in &mut m {
*x *= k;
}
}
for x in &mut m {
if !ffi::svd_element_ok(*x) {
*x = 0.0;
}
}
Inputs {
op: "svd".into(),
rows,
cols,
m: m.into_iter().map(num).collect(),
v: Vec::new(),
i: 0,
j: 0,
s: 0.0,
}
}
/// A 3x3 block of the kinds `SmBrep::Transform(double[4][4])` decomposes.
fn transform_block(rng: &mut Rng) -> Vec<f64> {
let (a, b) = (rng.uniform(-3.2, 3.2), rng.uniform(-3.2, 3.2));
let (ca, sa, cb, sb) = (a.cos(), a.sin(), b.cos(), b.sin());
// a rotation about Z then X
let rot = [ca, sa, 0.0, -sa * cb, ca * cb, sb, sa * sb, -ca * sb, cb];
let mut d = [1.0, 1.0, 1.0];
match rng.below(6) {
0 => {}
1 => d = [2.0, 2.0, 2.0],
2 => d = [rng.uniform(0.1, 5.0), rng.uniform(0.1, 5.0), rng.uniform(0.1, 5.0)],
3 => d = [-1.0, 1.0, 1.0],
4 => d = [1.0, 1.0, 10f64.powi(-(rng.range_i64(3, 15) as i32))],
_ => d = [1.0, 1.0, 0.0],
}
let mut m: Vec<f64> = rot
.chunks(3)
.zip(d)
.flat_map(|(row, dk)| row.iter().map(move |x| x * dk))
.collect();
if rng.chance(0.4) {
// a shear term
let k = rng.below(9) as usize;
if let Some(x) = m.get_mut(k) {
*x += rng.uniform(-1.0, 1.0);
}
}
m
}
impl UnitDiff for Runner {
fn id(&self) -> &'static str {
"geom.matrix"
@@ -487,7 +575,19 @@ impl UnitDiff for Runner {
seeds: ctx.seeds.max(MIN_GENERATED),
..ctx.clone()
};
default_cases(self, &ctx)
let mut cases = default_cases(self, &ctx)?;
for i in 0..SVD_GENERATED {
let seed = CaseContext::seed_for("geom.matrix.svd", i);
let mut rng = Rng::new(seed);
cases.push(Case {
unit: self.id().to_string(),
name: format!("svd{i:04}"),
source: crate::case::CaseSource::Generated { seed },
inputs: serde_json::to_value(gen_svd(&mut rng))
.map_err(|e| HarnessError::new(e.to_string()))?,
});
}
Ok(cases)
}
fn run_cpp(&self, case: &Case) -> Result<Output, HarnessError> {
@@ -8,13 +8,20 @@
//! edge curve, every edge use's UV trim curve and every face surface (type code,
//! natural domain, evaluated points). Inputs the native reader cannot build are
//! `rejected`, not failed.
//!
//! A case with a `matrix` runs `SmBrep::Transform(const double[4][4], bScaleTol)`
//! (`topology::brep::mirror::transform_matrix`) instead of the mirror, compared the same
//! way: the fixture file `transform_matrix.json` (identity, rotation, uniform and
//! non-uniform scale, shear, mirror, near-singular, singular, huge and tiny scale) and
//! [`MATRIX_GENERATED`] generated `mat####` cases, seeded apart from the `gen####`
//! mirror stream.
use serde::{Deserialize, Serialize};
use serde_json::Value;
use solidcore::geom::axis2_placement::Axis2Placement;
use solidcore::kernel::Kernel;
use solidcore::math::Vec3;
use solidcore::topology::brep::mirror::mirror;
use solidcore::topology::brep::mirror::{mirror, transform_matrix};
use solidcore::topology::brep::p3::notify_serial;
use solidcore::topology::Brep;
use solidcore::{Result as SolidResult, SolidError};
@@ -57,8 +64,16 @@ struct Inputs {
origin: [f64; 3],
x: [f64; 3],
y: [f64; 3],
/// When set, `SmBrep::Transform` by this 4x4 matrix (rows; the translation is row
/// 3) instead of the mirror.
matrix: Option<[[f64; 4]; 4]>,
/// `bScaleTol` of the matrix transform (upstream's default `TRUE`).
scale_tol: bool,
}
/// Generated matrix-transform cases per run.
pub const MATRIX_GENERATED: usize = 96;
impl Default for Inputs {
fn default() -> Self {
Self {
@@ -68,6 +83,8 @@ impl Default for Inputs {
origin: [0.; 3],
x: [1., 0., 0.],
y: [0., 1., 0.],
matrix: None,
scale_tol: true,
}
}
}
@@ -109,7 +126,19 @@ fn path(i: &Inputs) -> Result<String, HarnessError> {
/// The native output, or `Err(output)` when the FFI guard refuses the input.
fn native(i: &Inputs, op: u8) -> Result<Result<Vec<f64>, Output>, HarnessError> {
match ffi::eval(op, &path(i)?, &block(i), &args(i)) {
let r = match (op, i.matrix) {
(1, Some(m)) => {
let mut a = [0.0; 18];
a[0] = f64::from(i.editing);
a[1] = f64::from(i.scale_tol);
for (slot, x) in a.iter_mut().skip(2).zip(m.iter().flatten()) {
*slot = *x;
}
ffi::transform(&path(i)?, &block(i), &a)
}
_ => ffi::eval(op, &path(i)?, &block(i), &args(i)),
};
match r {
Ok(raw) if raw.first() == Some(&-999.) => Ok(Err(Output::rejected(format!(
"fixture constructor status {}",
raw.get(1).copied().unwrap_or(-1.)
@@ -251,7 +280,11 @@ impl UnitDiff for Runner {
b.set_editing_enabled(i.editing);
let plane = Axis2Placement::new(v(i.origin), v(i.x), v(i.y));
let mut run = || -> SolidResult<Vec<f64>> {
let st = status(mirror(&mut k, &mut b, &plane))?;
let r = match &i.matrix {
Some(m) => transform_matrix(&mut k, &mut b, m, i.scale_tol),
None => mirror(&mut k, &mut b, &plane),
};
let st = status(r)?;
let mut record = vec![st];
record.extend(geometry(&k, &b, true)?);
Ok(record)
@@ -319,10 +352,86 @@ impl UnitDiff for Runner {
fn cases(&self, ctx: &CaseContext) -> Result<Vec<Case>, HarnessError> {
let mut ctx = ctx.clone();
ctx.seeds = ctx.seeds.max(128);
crate::runner::default_cases(self, &ctx)
let mut cases = crate::runner::default_cases(self, &ctx)?;
for n in 0..MATRIX_GENERATED {
let seed = CaseContext::seed_for("topology.brep_mirror.transform_matrix", n);
let mut rng = Rng::new(seed);
let Some(inputs) = generate_matrix(&mut rng) else {
break;
};
cases.push(Case {
unit: self.id().to_string(),
name: format!("mat{n:04}"),
source: crate::case::CaseSource::Generated { seed },
inputs,
});
}
Ok(cases)
}
}
/// A matrix-transform case: a generated 3x3 block (rotation times a scale of one of
/// the kinds below, sometimes sheared) and translation, on one of [`FILES`].
fn generate_matrix(rng: &mut Rng) -> Option<Value> {
let file = FILES.get(rng.below(FILES.len() as u64) as usize)?;
let (a, c) = (rng.uniform(-3.2, 3.2), rng.uniform(-3.2, 3.2));
let (ca, sa, cc, sc) = (a.cos(), a.sin(), c.cos(), c.sin());
// a rotation about Z then X
let rot = [[ca, sa, 0.], [-sa * cc, ca * cc, sc], [sa * sc, -ca * sc, cc]];
let d = match rng.below(9) {
0 => [1., 1., 1.],
1 => {
let s = rng.uniform(0.2, 5.0);
[s, s, s]
}
2 | 3 => [
rng.uniform(0.2, 5.0),
rng.uniform(0.2, 5.0),
rng.uniform(0.2, 5.0),
],
4 => [-1., 1., 1.],
5 => [-rng.uniform(0.2, 5.0), rng.uniform(0.2, 5.0), 1.],
6 => [1., 1., 10f64.powi(-(rng.range_i64(3, 10) as i32))],
7 => {
let s = 10f64.powi(rng.range_i64(-6, 6) as i32);
[s, s, s]
}
_ => [1., 1., 0.],
};
let mut m = [[0.; 4]; 4];
for ((out, row), dr) in m.iter_mut().zip(rot).zip(d) {
for (o, x) in out.iter_mut().zip(row) {
*o = x * dr;
}
}
if rng.chance(0.4) {
// a shear term
let (r, k) = (rng.below(3) as usize, rng.below(3) as usize);
let shear = rng.uniform(-1.0, 1.0);
if let Some(x) = m.get_mut(r).and_then(|row| row.get_mut(k)) {
*x += shear;
}
}
let scale = 10f64.powi(rng.range_i64(-1, 2) as i32);
let t = [
rng.uniform(-10., 10.) * scale,
rng.uniform(-10., 10.) * scale,
rng.uniform(-10., 10.) * scale,
1.,
];
if let Some(row) = m.get_mut(3) {
*row = t;
}
let i = Inputs {
file: (*file).into(),
matrix: Some(m),
scale_tol: rng.chance(0.8),
editing: rng.chance(0.9),
..Inputs::default()
};
serde_json::to_value(i).ok()
}
/// `Face::shell` and `Shell::faces` have no upstream counterpart: the import derives them
/// from the first face use in walk order, the port keeps them with the face's first
/// face use (`swap_reflection_shells`). Both are recomputed here from the face's own
@@ -66,6 +66,22 @@ pub enum Op {
/// `Compute3x3EigenVectors`: values = 3 eigenvalues, 3 vectors (xyz each),
/// matrix after. `SM_ERR` for a non-symmetric matrix (upstream output undefined).
Eigen3 = 10,
/// `SVDecompose(W, V)`: values = `W` (`cols`), `V` (`cols x cols`, row-major),
/// matrix after (`U`). Every element must be zero or of magnitude in
/// [`SVD_MIN_ABS`, `SVD_MAX_ABS`]: a NaN or infinite norm makes upstream's
/// splitting test read `rW[-1]`.
Svd = 11,
}
/// Smallest non-zero element magnitude [`Op::Svd`] accepts.
pub const SVD_MIN_ABS: f64 = 1e-100;
/// Largest element magnitude [`Op::Svd`] accepts.
pub const SVD_MAX_ABS: f64 = 1e100;
/// True when `x` is zero or a finite magnitude in [`SVD_MIN_ABS`, `SVD_MAX_ABS`].
#[must_use]
pub fn svd_element_ok(x: f64) -> bool {
x == 0.0 || (SVD_MIN_ABS..=SVD_MAX_ABS).contains(&x.abs())
}
/// Output of [`op`].
@@ -106,6 +122,7 @@ pub fn op(
Op::Resize => (1..=MAX_DIM as i64).contains(&i) && (1..=MAX_DIM as i64).contains(&j),
Op::SolveLu => rows == cols && v.len() == rows,
Op::Eigen3 => rows == 3 && cols == 3,
Op::Svd => m.iter().all(|&x| svd_element_ok(x)),
_ => true,
};
if !ok {
@@ -6,7 +6,8 @@
//! `Mirror` and returns `record_len record* 0 graph`, the record being the status, then
//! samples of every edge curve, its edge uses' UV trim curves and every face surface
//! (type, natural interval or domain, five or nine evaluated points) and the raw graph
//! of the mirrored B-rep. A failed fixture build returns `-999 status`.
//! of the mirrored B-rep. A failed fixture build returns `-999 status`. Op 3 ([`transform`])
//! runs `SmBrep::Transform(const double[4][4], SmBoolean)` instead, with the same output.
use super::topology_edgeuse_p1::{brep_block, Reader};
use crate::{call_with_buffer, FfiError};
@@ -40,6 +41,32 @@ pub fn eval(op: u8, path: &str, block: &[f64], args: &[f64; 10]) -> Result<Vec<f
if op > 2 || !args.iter().all(|x| x.is_finite()) || !matches!(args[0], 0.0 | 1.0) {
return Err(bad("mirror placement or editing flag"));
}
run(op, path, block, args)
}
/// Op 3: `SmBrep::Transform(const double[4][4], bScaleTol)` on the fixture, output as
/// op 1 of [`eval`].
///
/// `args`: editing flag, `bScaleTol` (each 0 or 1), then the 4x4 matrix row-major (the
/// 3x3 block is rows 0..3 / columns 0..3, the translation row 3). Every 3x3 element must
/// be zero or of magnitude in `[1e-100, 1e100]`, as [`super::geom_matrix::Op::Svd`]
/// requires (a NaN or infinite norm makes `SVDecompose` read `rW[-1]`), and the
/// translation finite.
pub fn transform(path: &str, block: &[f64], args: &[f64; 18]) -> Result<Vec<f64>, FfiError> {
let flag = |x: f64| matches!(x, 0.0 | 1.0);
let block3 = args
.get(2..14)
.is_some_and(|m| m.chunks(4).all(|row| row.iter().take(3).all(|&x| super::geom_matrix::svd_element_ok(x))));
if !args.iter().all(|x| x.is_finite()) || !flag(args[0]) || !flag(args[1]) || !block3 {
return Err(FfiError::InvalidArgument(
"transform matrix or editing / scale-tolerance flag".into(),
));
}
run(3, path, block, args)
}
fn run(op: u8, path: &str, block: &[f64], args: &[f64]) -> Result<Vec<f64>, FfiError> {
let bad = |m: &str| FfiError::InvalidArgument(m.into());
let mut rd = Reader {
input: block,
at: 0,
@@ -57,7 +84,7 @@ pub fn eval(op: u8, path: &str, block: &[f64], args: &[f64; 10]) -> Result<Vec<f
block.as_ptr(),
block.len() as i64,
args.as_ptr(),
10,
args.len() as i64,
out.as_mut_ptr(),
out.len() as i64,
&raw mut *len,
+416 -1
View File
@@ -25,9 +25,12 @@
//! | [`Matrix::solve_linear_system`] | `SolveLinearSystem` |
//! | [`Matrix::eigen_system_jacobi`] | `EigenSystemJacobiSolver` |
//! | [`Matrix::compute_3x3_eigen_vectors`] | `Compute3x3EigenVectors` |
//! | [`Matrix::transpose`] | `Transpose` |
//! | [`Matrix::sv_decompose`], [`SVD_MAX_ITERATIONS`] | `SVDecompose`, its `lMaxIters` |
//! | `pythag` (private) | `smos_Pythag` (`SmMath.h`) |
//! | `sort_eigen_system` (private) | `sm_SortEigenSystem` |
//!
//! The other `SmMatrix` members (`Add`, `Multiply`, `Invert`, `Transpose`, SVD,
//! The other `SmMatrix` members (`Add`, `Multiply`, `Invert`, `SolveWithSVDMatrices`,
//! `Rank`, `MakeRowEchelon`, `SwapRows`, 3x3 determinant / eigenvalues, ...) are not
//! reached by any upstream test entry point and are not ported (see
//! `harness/callgraph/units.json`).
@@ -64,6 +67,10 @@
//! there, and `InvalidInput` for a matrix that is not 3x3 (upstream reads out of
//! bounds). When the solver does not converge, the unconverged values are
//! returned with success, as upstream.
//! - `SVDecompose` returns `W` and `V` instead of filling out-parameters, and omits
//! its `SM_DEBUG_CODE` self-check. Where a NaN norm makes its splitting test read
//! `rW[-1]` (undefined upstream) it returns `InvalidInput`. Its QR loop is bounded
//! by `lMaxIters` (30) per singular value, as upstream; no recursion.
use crate::error::{Result, SolidError};
use crate::geom::tol::is_zero;
@@ -94,6 +101,30 @@ fn sm_err(operation: &'static str) -> SolidError {
}
}
/// `lMaxIters` of `SVDecompose`: QR sweeps per singular value.
pub const SVD_MAX_ITERATIONS: usize = 30;
/// `smos_Pythag(a, b)`: `sqrt(a*a + b*b)` without destructive overflow or underflow,
/// signed like its larger operand (the macro of `SmMath.h`).
fn pythag(a: f64, b: f64) -> f64 {
if a.abs() > b.abs() {
a * sqrt(1.0 + (b / a) * (b / a))
} else if b != 0.0 {
b * sqrt(1.0 + (a / b) * (a / b))
} else {
0.0
}
}
fn vget(v: &[f64], i: usize, op: &'static str) -> Result<f64> {
v.get(i).copied().ok_or_else(|| out_of_range(op))
}
fn vset(v: &mut [f64], i: usize, x: f64, op: &'static str) -> Result<()> {
*v.get_mut(i).ok_or_else(|| out_of_range(op))? = x;
Ok(())
}
/// A dense `rows x cols` matrix stored row-major (`SmMatrix`).
#[derive(Debug, Clone, PartialEq)]
pub struct Matrix {
@@ -686,6 +717,390 @@ impl Matrix {
};
Ok((values, [col(0)?, col(1)?, col(2)?]))
}
/// Transposes the matrix in place (`Transpose`): a square matrix swaps its
/// off-diagonal pairs, any other becomes the `cols x rows` transpose.
pub fn transpose(&mut self) -> Result<()> {
const OP: &str = "geom::matrix::Matrix::transpose";
let (m, n) = (self.rows, self.cols);
if m == n {
for ii in 0..m {
for jj in ii + 1..m {
let temp = self.get(ii, jj, OP)?;
let other = self.get(jj, ii, OP)?;
self.put(ii, jj, other, OP)?;
self.put(jj, ii, temp, OP)?;
}
}
} else {
let mut t = Self::try_new(n, m)?;
for ii in 0..m {
for jj in 0..n {
t.put(jj, ii, self.get(ii, jj, OP)?, OP)?;
}
}
*self = t;
}
Ok(())
}
/// Singular value decomposition `A = U W Vᵀ` of this `m x n` matrix (`m >= n`;
/// `SVDecompose`, after "Numerical Recipes in C"): the matrix is overwritten by
/// `U` (`m x n`); returns the diagonal of `W` (`n` values) and `V` (`n x n`, not
/// `Vᵀ`).
///
/// The QR sweep runs at most [`SVD_MAX_ITERATIONS`] times per singular value; when a
/// value has not converged by then the decomposition carries on with the next one
/// and returns success, as upstream (its `its == lMaxIters` error test is never
/// reached by the zero-based loop).
///
/// # Errors
///
/// `SM_ERR` when `m < n` or the matrix is empty. `InvalidInput` where upstream
/// reads outside its arrays (a NaN norm makes the splitting test fall through to
/// `rv1[-1]`).
pub fn sv_decompose(&mut self) -> Result<(Vec<f64>, Self)> {
const OP: &str = "geom::matrix::Matrix::sv_decompose";
let m = self.rows;
let n = self.cols;
// check input: m >= n
if m < n || m == 0 || n == 0 {
return Err(sm_err(OP));
}
let mut w = vec![0.0; n];
let mut v = Self::try_new(n, n)?;
let mut rv1 = vec![0.0; n];
let a = self;
let mut anorm: f64 = 0.0;
let mut g: f64 = 0.0;
let mut scale: f64 = 0.0;
let (mut c, mut f, mut h, mut s);
let (mut x, mut y, mut z);
// Householder reduction to bidiagonal form
let mut ll;
for i1 in 0..n {
ll = i1 + 1;
vset(&mut rv1, i1, scale * g, OP)?;
g = 0.0;
s = 0.0;
scale = 0.0;
if i1 < m {
for kk in i1..m {
scale += a.get(kk, i1, OP)?.abs();
}
if scale != 0.0 {
for kk in i1..m {
let q = a.get(kk, i1, OP)? / scale;
a.put(kk, i1, q, OP)?;
s += q * q;
}
f = a.get(i1, i1, OP)?;
g = -(if f >= 0.0 { sqrt(s) } else { -sqrt(s) });
h = f * g - s;
a.put(i1, i1, f - g, OP)?;
if i1 != n - 1 {
for j1 in ll..n {
s = 0.0;
for kk in i1..m {
s += a.get(kk, i1, OP)? * a.get(kk, j1, OP)?;
}
f = s / h;
for kk in i1..m {
let q = a.get(kk, j1, OP)? + f * a.get(kk, i1, OP)?;
a.put(kk, j1, q, OP)?;
}
}
}
for kk in i1..m {
let q = a.get(kk, i1, OP)? * scale;
a.put(kk, i1, q, OP)?;
}
}
}
vset(&mut w, i1, scale * g, OP)?;
g = 0.0;
s = 0.0;
scale = 0.0;
if i1 < m && i1 != n - 1 {
for kk in ll..n {
scale += a.get(i1, kk, OP)?.abs();
}
if scale != 0.0 {
for kk in ll..n {
let q = a.get(i1, kk, OP)? / scale;
a.put(i1, kk, q, OP)?;
s += q * q;
}
f = a.get(i1, ll, OP)?;
g = -(if f >= 0.0 { sqrt(s) } else { -sqrt(s) });
h = f * g - s;
a.put(i1, ll, f - g, OP)?;
for kk in ll..n {
vset(&mut rv1, kk, a.get(i1, kk, OP)? / h, OP)?;
}
if i1 != m - 1 {
for j1 in ll..m {
s = 0.0;
for kk in ll..n {
s += a.get(j1, kk, OP)? * a.get(i1, kk, OP)?;
}
for kk in ll..n {
let q = a.get(j1, kk, OP)? + s * vget(&rv1, kk, OP)?;
a.put(j1, kk, q, OP)?;
}
}
}
for kk in ll..n {
let q = a.get(i1, kk, OP)? * scale;
a.put(i1, kk, q, OP)?;
}
}
}
let norm = vget(&w, i1, OP)?.abs() + vget(&rv1, i1, OP)?.abs();
anorm = if anorm > norm { anorm } else { norm };
}
// Accumulation of right hand transformations
ll = 0;
for i1 in (0..n).rev() {
if i1 < n - 1 {
if g != 0.0 {
for j1 in ll..n {
let q = (a.get(i1, j1, OP)? / a.get(i1, ll, OP)?) / g;
v.put(j1, i1, q, OP)?;
}
for j1 in ll..n {
s = 0.0;
for kk in ll..n {
s += a.get(i1, kk, OP)? * v.get(kk, j1, OP)?;
}
for kk in ll..n {
let q = v.get(kk, j1, OP)? + s * v.get(kk, i1, OP)?;
v.put(kk, j1, q, OP)?;
}
}
}
for j1 in ll..n {
v.put(j1, i1, 0.0, OP)?;
v.put(i1, j1, 0.0, OP)?;
}
}
v.put(i1, i1, 1.0, OP)?;
g = vget(&rv1, i1, OP)?;
ll = i1;
}
// Accumulation of left hand transformations
for i1 in (0..n).rev() {
ll = i1 + 1;
g = vget(&w, i1, OP)?;
if i1 < n - 1 {
for j1 in ll..n {
a.put(i1, j1, 0.0, OP)?;
}
}
if g == 0.0 {
for j1 in i1..m {
a.put(j1, i1, 0.0, OP)?;
}
} else {
g = 1.0 / g;
if i1 != n - 1 {
for j1 in ll..n {
s = 0.0;
for kk in ll..m {
s += a.get(kk, i1, OP)? * a.get(kk, j1, OP)?;
}
f = (s / a.get(i1, i1, OP)?) * g;
for kk in i1..m {
let q = a.get(kk, j1, OP)? + f * a.get(kk, i1, OP)?;
a.put(kk, j1, q, OP)?;
}
}
}
for j1 in i1..m {
let q = a.get(j1, i1, OP)? * g;
a.put(j1, i1, q, OP)?;
}
}
let q = a.get(i1, i1, OP)? + 1.0;
a.put(i1, i1, q, OP)?;
}
// Diagonalization of the bidiagonal form
let mut nm: usize;
for kk in (0..n).rev() {
// loop over allowed iterations
for its in 0..SVD_MAX_ITERATIONS {
let mut flag = true;
// test for splitting; rv1[0] is always 0
let mut l_int = kk;
loop {
ll = l_int;
nm = ll.wrapping_sub(1);
if vget(&rv1, ll, OP)?.abs() + anorm == anorm {
flag = false;
break;
}
if vget(&w, nm, OP)?.abs() + anorm == anorm {
break;
}
if l_int == 0 {
break;
}
l_int -= 1;
}
if flag {
// cancellation of rv1[ll] if ll > 1 (upstream also sets c = 0, which
// is never read)
s = 1.0;
for i1 in ll..=kk {
f = s * vget(&rv1, i1, OP)?;
if f.abs() + anorm != anorm {
g = vget(&w, i1, OP)?;
h = pythag(f, g);
vset(&mut w, i1, h, OP)?;
h = 1.0 / h;
c = g * h;
s = -f * h;
for j1 in 0..m {
y = a.get(j1, nm, OP)?;
z = a.get(j1, i1, OP)?;
a.put(j1, nm, y * c + z * s, OP)?;
a.put(j1, i1, z * c - y * s, OP)?;
}
}
}
}
z = vget(&w, kk, OP)?;
if ll == kk {
// Convergence.
if z < 0.0 {
// Singular value is made nonnegative
vset(&mut w, kk, -z, OP)?;
for j1 in 0..n {
let q = -v.get(j1, kk, OP)?;
v.put(j1, kk, q, OP)?;
}
}
break;
}
// upstream's `its == lMaxIters` error exit: unreachable in the
// zero-based loop, kept for the same control flow
if its == SVD_MAX_ITERATIONS {
return Err(sm_err(OP));
}
// shift from bottom 2-by-2 minor
x = vget(&w, ll, OP)?;
nm = kk.wrapping_sub(1);
y = vget(&w, nm, OP)?;
g = vget(&rv1, nm, OP)?;
h = vget(&rv1, kk, OP)?;
f = ((y - z) * (y + z) + (g - h) * (g + h)) / (2.0 * h * y);
g = (f * f + 1.0).sqrt();
f = ((x - z) * (x + z)
+ h * ((y / (f + if f >= 0.0 { g.abs() } else { -g.abs() })) - h))
/ x;
// next QR transformation
c = 1.0;
s = 1.0;
for j1 in ll..=nm {
let i1 = j1 + 1;
g = vget(&rv1, i1, OP)?;
y = vget(&w, i1, OP)?;
h = s * g;
g *= c;
z = pythag(f, h);
vset(&mut rv1, j1, z, OP)?;
c = f / z;
s = h / z;
f = x * c + g * s;
g = g * c - x * s;
h = y * s;
y *= c;
for jj in 0..n {
x = v.get(jj, j1, OP)?;
z = v.get(jj, i1, OP)?;
v.put(jj, j1, x * c + z * s, OP)?;
v.put(jj, i1, z * c - x * s, OP)?;
}
z = pythag(f, h);
// rotation can be arbitrary if z = 0
vset(&mut w, j1, z, OP)?;
if z != 0.0 {
z = 1.0 / z;
c = f * z;
s = h * z;
}
f = (c * g) + (s * y);
x = (c * y) - (s * g);
for jj in 0..m {
y = a.get(jj, j1, OP)?;
z = a.get(jj, i1, OP)?;
a.put(jj, j1, y * c + z * s, OP)?;
a.put(jj, i1, z * c - y * s, OP)?;
}
}
vset(&mut rv1, ll, 0.0, OP)?;
vset(&mut rv1, kk, f, OP)?;
vset(&mut w, kk, x, OP)?;
}
}
Ok((w, v))
}
}
/// Sorts eigenvalues ascending, swapping the matching eigenvector columns
+240 -2
View File
@@ -6,6 +6,8 @@
//! | Rust | Upstream |
//! |-------------------------------------------------|--------------------------------------------|
//! | [`mirror`], [`mirror_with_notifications`] | `SmBrep::Mirror` |
//! | [`transform_matrix`], [`transform_matrix_with_notifications`] | `SmBrep::Transform(const double[4][4], SmBoolean)` |
//! | `fold_status` (private) | `eStat = std::max(..., eStat)` |
//! | `Axis2Placement::mirror_point` | `SmAxis2Placement::MirrorPoint` |
//! | `Axis2Placement::mirror_vector` | `SmAxis2Placement::MirrorVector` |
//! | `Curve::create_mirror_curve` | `SmCurve::CreateMirrorCurve` (virtual) |
@@ -30,10 +32,19 @@
//! alone). The NURBs of analytic surfaces without one are made before the infinite region
//! search (`make_missing_nurbs`), where upstream's bounding box makes them lazily. The
//! debug drawing code is omitted. No recursion is introduced here.
//!
//! [`transform_matrix`] runs upstream's steps and branches in order: the low-work
//! test, `SVDecompose` (status ignored upstream; a failure here is `SM_ERR` without
//! touching the B-rep, where upstream reads the arrays the failed call left), then
//! the `Transform` / `Mirror` sequence, every status folded with `std::max` (the error
//! with the largest code is returned; one without an exact status counts as `SM_ERR`).
//! A failed `AngleBetween` reads as angle `0`, as upstream sets it.
use crate::error::{ignore_status, Result, SolidError};
use crate::geom::axis2_placement::Axis2Placement;
use crate::geom::matrix::Matrix;
use crate::kernel::Kernel;
use crate::math::{Vec3, EFF_ZERO, EFF_ZERO_RAD, HALF_PI};
use crate::query::brep::BrepQuery;
use crate::surface::surface::p2::SurfaceKind;
use crate::topology::brep::p1::{
@@ -44,7 +55,7 @@ use crate::topology::brep::p5::{
brep_faces, brep_regions, brep_vertices, find_and_set_infinite_region, notify_edge_pre_edit,
notify_face_pre_edit, notify_pre_edit, notify_vertex_pre_edit, plain_edit, replace_surface,
set_edge_curve_owner, set_joined_edge_curve, set_vertex_point, transform_owner_key,
with_notifications, EditNotifyEnv,
transform_with_notifications, with_notifications, EditNotifyEnv,
};
use crate::topology::Brep;
use crate::EntityRef;
@@ -252,10 +263,237 @@ fn make_missing_nurbs(k: &Kernel, b: &mut Brep) -> Result<()> {
Ok(())
}
/// Keeps the status of a step whose `SmStatus` upstream folds in with `std::max`:
/// `worst` holds the error with the largest status code so far (the first one on a
/// tie). Unported dependencies and preparation requests still propagate.
fn fold_status(worst: &mut Option<SolidError>, r: Result<()>) -> Result<()> {
let Err(e) = r else { return Ok(()) };
if e.is_unported() || matches!(e, SolidError::PreparationRequired { .. }) {
return Err(e);
}
let code = |e: &SolidError| e.sm_status().unwrap_or(crate::error::SM_ERR);
if worst.as_ref().is_none_or(|w| code(&e) > code(w)) {
*worst = Some(e);
}
Ok(())
}
/// `SmBrep::Transform(const double[4][4], bScaleTol)`: applies the affine map
/// `p' = p·A + t` to a B-rep, where `A` is the upper-left 3x3 block of `matrix`
/// (rows `matrix[0..3]`, columns `0..3`) and `t` its fourth row (row-vector
/// convention). The fourth column is ignored.
///
/// Without shear or reflection (rows 0 and 1 perpendicular within `SM_EFF_ZERO` and
/// row 0 × row 1 along row 2 within `SM_EFF_ZERO_RAD`) this is one
/// [`transform`](crate::topology::brep::p5::transform) by the unit rows scaled by
/// their lengths. Otherwise `A = V D Uᵀ` by [`Matrix::sv_decompose`] and the B-rep is
/// transformed by `V` (scaled by its row lengths), mirrored when `V` is left-handed,
/// then transformed by `Uᵀ` scaled by `D`, mirrored and translated when `Uᵀ` is
/// left-handed. A non-uniform scale turns every geometry into NURBS. `scale_tolerances`
/// is passed to every scaling transform.
///
/// # Errors
///
/// As upstream, every step runs even when an earlier one fails, and the error with
/// the largest upstream status is returned; those of
/// [`transform`](crate::topology::brep::p5::transform), [`mirror`] and
/// [`Matrix::sv_decompose`] (whose status upstream ignores).
pub fn transform_matrix(
k: &mut Kernel,
b: &mut Brep,
matrix: &[[f64; 4]; 4],
scale_tolerances: bool,
) -> Result<()> {
transform_matrix_with_notifications(k, b, matrix, scale_tolerances, None)
}
/// [`transform_matrix`] notifying every step through the caller's environment
/// `notifications` (each step is its own plain edit under the kernel's identities when
/// `None`, as upstream's separate `Transform` and `Mirror` calls).
///
/// # Errors
///
/// Those of [`transform_matrix`].
pub fn transform_matrix_with_notifications(
k: &mut Kernel,
b: &mut Brep,
matrix: &[[f64; 4]; 4],
scale_tolerances: bool,
mut notifications: Option<&mut EditNotifyEnv<'_>>,
) -> Result<()> {
const OP: &str = "topology::brep::mirror::transform_matrix";
let [m0, m1, m2, m3] = matrix;
// Setup 3x3 for SVD decomposition
let row0 = Vec3::new(m0[0], m0[1], m0[2]);
let row1 = Vec3::new(m1[0], m1[1], m1[2]);
let row2 = Vec3::new(m2[0], m2[1], m2[2]);
// Set up displacement
let disp = Vec3::new(m3[0], m3[1], m3[2]);
// Low work case: No sheer or reflection?
if row0.dot(row1).abs() < EFF_ZERO {
// a failed AngleBetween reports 0
let angle_check = row0.cross(row1).angle_between(row2).unwrap_or(0.0);
if angle_check <= EFF_ZERO_RAD {
let row0_len = row0.length();
let row1_len = row1.length();
let row2_len = row2.length();
let input_scale = Vec3::new(row0_len, row1_len, row2_len);
let input_transform = Axis2Placement::new(disp, row0 / row0_len, row1 / row1_len);
// Use SmBrep::Transform with no shear/reflection
return transform_with_notifications(
k,
b,
&input_transform,
Some(input_scale),
scale_tolerances,
notifications,
);
}
}
let mut sv = Matrix::from_row_major(
3,
3,
vec![
row0.x, row0.y, row0.z, row1.x, row1.y, row1.z, row2.x, row2.y, row2.z,
],
)?;
// Decompose matrix (upstream ignores the status; `sU` and `sDDiag` then hold
// whatever the failed call left)
let (d_diag, mut ut) = match sv.sv_decompose() {
Ok(r) => r,
Err(e) if e.is_unported() => return Err(e),
Err(_) => return Err(SolidError::sm_err(OP)),
};
ut.transpose()?;
let row = |m: &Matrix, r: usize| -> Result<Vec3> {
let v = m
.row(r)
.ok_or_else(|| SolidError::port_invariant(OP, "3x3 matrix row"))?;
match v {
[x, y, z] => Ok(Vec3::new(*x, *y, *z)),
_ => Err(SolidError::port_invariant(OP, "3x3 matrix row")),
}
};
let d = |i: usize| -> Result<f64> {
d_diag
.get(i)
.copied()
.ok_or_else(|| SolidError::port_invariant(OP, "3 singular values"))
};
// Build Axis2Placement and scaling vectors for transforms
let vx = row(&sv, 0)?;
let vy = row(&sv, 1)?;
let vz = row(&sv, 2)?;
let vx_len = vx.length();
let vy_len = vy.length();
let vz_len = vz.length();
let mut transform_v = Axis2Placement::IDENTITY;
let scale_v = Vec3::new(vx_len, vy_len, vz_len);
let origin = Vec3::new(0.0, 0.0, 0.0);
// upstream ignores SetCanonical's status
let _ = transform_v.set_canonical(origin, vx / vx_len, vy / vy_len);
let utx = row(&ut, 0)?;
let uty = row(&ut, 1)?;
let utz = row(&ut, 2)?;
let utx_len = utx.length();
let uty_len = uty.length();
let utz_len = utz.length();
let mut transform_ut = Axis2Placement::IDENTITY;
let scale_ut = Vec3::new(utx_len * d(0)?, uty_len * d(1)?, utz_len * d(2)?);
// Transforms below are rigid, need to check if we require a mirror to fix the
// handedness of the orthonormal decomposition.
let angle_v = vx.cross(vy).angle_between(vz).unwrap_or(0.0);
let mirror_v = angle_v > HALF_PI;
let angle_ut = utx.cross(uty).angle_between(utz).unwrap_or(0.0);
let mirror_ut = angle_ut > HALF_PI;
let mut worst = None;
// Transform the brep, applying the mirrors as necessary. Note that the translation
// comes at the very end.
fold_status(
&mut worst,
transform_with_notifications(
k,
b,
&transform_v,
Some(scale_v),
scale_tolerances,
notifications.as_deref_mut(),
),
)?;
if mirror_v {
fold_status(
&mut worst,
mirror_with_notifications(k, b, &transform_v, notifications.as_deref_mut()),
)?;
}
if mirror_ut {
let _ = transform_ut.set_canonical(origin, utx / utx_len, uty / uty_len);
fold_status(
&mut worst,
transform_with_notifications(
k,
b,
&transform_ut,
Some(scale_ut),
scale_tolerances,
notifications.as_deref_mut(),
),
)?;
fold_status(
&mut worst,
mirror_with_notifications(k, b, &transform_ut, notifications.as_deref_mut()),
)?;
let v_x = Vec3::new(1.0, 0.0, 0.0);
let v_y = Vec3::new(0.0, 1.0, 0.0);
let _ = transform_ut.set_canonical(disp, v_x, v_y);
// the defaults: no scale, bScaleTol = TRUE
fold_status(
&mut worst,
transform_with_notifications(k, b, &transform_ut, None, true, notifications),
)?;
} else {
// transform without mirror
let _ = transform_ut.set_canonical(disp, utx / utx_len, uty / uty_len);
fold_status(
&mut worst,
transform_with_notifications(
k,
b,
&transform_ut,
Some(scale_ut),
scale_tolerances,
notifications,
),
)?;
}
// all done
worst.map_or(Ok(()), Err)
}
#[cfg(test)]
mod tests {
use super::*;
use crate::math::Vec3;
use crate::primitive::api_primitives::create_box;
use crate::topology::brep::p1;
@@ -0,0 +1,123 @@
//! Issue #1570: `SmBrep::Transform(const double[4][4], bScaleTol)` (and the
//! `SmMatrix::SVDecompose` it relies on) applies an arbitrary affine map `p' = p·A + t`
//! (row-vector convention) to a B-rep: shear, non-uniform scale and reflection included.
//!
//! The issue's cases on the unit cube `[0,1]^3`: every vertex maps to `v·A + t`, the
//! volume becomes `|det A|` and the shell stays outward (positive volume) under a
//! reflection.
#![allow(
clippy::unwrap_used,
clippy::expect_used,
clippy::panic,
clippy::indexing_slicing
)] // test helpers
use solidcore::geom::attribute::{AttributeStore, ObjectKey};
use solidcore::geom::matrix::Matrix;
use solidcore::mass::integrator::properties::{brep_volume, PropertyAttributes};
use solidcore::math::{Point3, Vec3};
use solidcore::primitive::api_primitives::create_box;
use solidcore::topology::brep::mirror::transform_matrix;
use solidcore::topology::brep::p1;
use solidcore::{Brep, Kernel};
fn volume(k: &Kernel, b: &mut Brep) -> f64 {
let store = AttributeStore::default();
let attrs = PropertyAttributes {
store: &store,
key_of: &|_| ObjectKey(0),
};
brep_volume(k, b, 1e-6, &attrs).unwrap()
}
fn apply(m: &[[f64; 4]; 4], p: Point3) -> Point3 {
Vec3::new(
p.x * m[0][0] + p.y * m[1][0] + p.z * m[2][0] + m[3][0],
p.x * m[0][1] + p.y * m[1][1] + p.z * m[2][1] + m[3][1],
p.x * m[0][2] + p.y * m[1][2] + p.z * m[2][2] + m[3][2],
)
}
/// Transforms the unit cube by `m` and checks every vertex and the volume.
fn check(m: [[f64; 4]; 4], expected_volume: f64) {
let mut k = Kernel::default();
let mut b = create_box(&mut k, Point3::new(0.0, 0.0, 0.0), 1.0, 1.0, 1.0).unwrap();
let vertices = p1::vertices(&b, None).unwrap();
let before: Vec<Point3> = vertices
.iter()
.map(|&v| b.vertex(v).unwrap().point)
.collect();
transform_matrix(&mut k, &mut b, &m, true).unwrap();
for (&v, &p) in vertices.iter().zip(&before) {
let got = b.vertex(v).unwrap().point;
let want = apply(&m, p);
assert!(
(got - want).length() < 1e-9,
"vertex {p:?}: got {got:?}, want {want:?}"
);
}
let vol = volume(&k, &mut b);
assert!(
(vol - expected_volume).abs() < 1e-6,
"volume {vol}, want {expected_volume}"
);
}
#[test]
fn issue_1570_shear_keeps_the_volume() {
let m = [
[1.0, 0.5, 0.0, 0.0],
[0.0, 1.0, 0.0, 0.0],
[0.0, 0.0, 1.0, 0.0],
[2.0, 3.0, 4.0, 1.0],
];
check(m, 1.0);
// the issue's spot check: (1,1,1) -> (3, 4.5, 5)
assert_eq!(
apply(&m, Point3::new(1.0, 1.0, 1.0)),
Point3::new(3.0, 4.5, 5.0)
);
}
#[test]
fn issue_1570_non_uniform_scale() {
check(
[
[2.0, 0.0, 0.0, 0.0],
[0.0, 1.0, 0.0, 0.0],
[0.0, 0.0, 1.0, 0.0],
[1.0, -2.0, 0.5, 1.0],
],
2.0,
);
}
#[test]
fn issue_1570_reflection_keeps_the_shell_outward() {
check(
[
[-1.0, 0.0, 0.0, 0.0],
[0.0, 1.0, 0.0, 0.0],
[0.0, 0.0, 1.0, 0.0],
[0.0, 0.0, 0.0, 1.0],
],
1.0,
);
}
#[test]
fn issue_1570_sv_decompose_reconstructs_the_matrix() {
let a = [1.0, 0.5, 0.0, 0.0, 1.0, 0.0, 0.3, -0.2, 2.0];
let mut u = Matrix::from_row_major(3, 3, a.to_vec()).unwrap();
let (w, v) = u.sv_decompose().unwrap();
for r in 0..3 {
for c in 0..3 {
let mut x = 0.0;
for (i, wi) in w.iter().enumerate() {
x += u.at(r, i).unwrap() * wi * v.at(c, i).unwrap();
}
assert!((x - a[r * 3 + c]).abs() < 1e-12, "A[{r}][{c}] = {x}");
}
}
assert!(w.iter().all(|&s| s >= 0.0));
}
+272
View File
@@ -0,0 +1,272 @@
{
"cases": [
{
"name": "svd_identity_3x3",
"upstream_ref": "source/SMLib/src/SmMatrix.cpp SmMatrix::SVDecompose",
"inputs": {
"op": "svd",
"rows": 3,
"cols": 3,
"m": [
1,
0,
0,
0,
1,
0,
0,
0,
1
]
}
},
{
"name": "svd_issue_1570_shear",
"upstream_ref": "source/SMLib/src/SmMatrix.cpp SmMatrix::SVDecompose",
"inputs": {
"op": "svd",
"rows": 3,
"cols": 3,
"m": [
1,
0.5,
0,
0,
1,
0,
0,
0,
1
]
}
},
{
"name": "svd_non_uniform_scale",
"upstream_ref": "source/SMLib/src/SmMatrix.cpp SmMatrix::SVDecompose",
"inputs": {
"op": "svd",
"rows": 3,
"cols": 3,
"m": [
2,
0,
0,
0,
1,
0,
0,
0,
1
]
}
},
{
"name": "svd_mirror_x",
"upstream_ref": "source/SMLib/src/SmMatrix.cpp SmMatrix::SVDecompose",
"inputs": {
"op": "svd",
"rows": 3,
"cols": 3,
"m": [
-1,
0,
0,
0,
1,
0,
0,
0,
1
]
}
},
{
"name": "svd_singular_z",
"upstream_ref": "source/SMLib/src/SmMatrix.cpp SmMatrix::SVDecompose",
"inputs": {
"op": "svd",
"rows": 3,
"cols": 3,
"m": [
1,
0,
0,
0,
1,
0,
0,
0,
0
]
}
},
{
"name": "svd_zero_3x3",
"upstream_ref": "source/SMLib/src/SmMatrix.cpp SmMatrix::SVDecompose",
"inputs": {
"op": "svd",
"rows": 3,
"cols": 3,
"m": [
0,
0,
0,
0,
0,
0,
0,
0,
0
]
}
},
{
"name": "svd_near_singular",
"upstream_ref": "source/SMLib/src/SmMatrix.cpp SmMatrix::SVDecompose",
"inputs": {
"op": "svd",
"rows": 3,
"cols": 3,
"m": [
1,
2,
3,
4,
5,
6,
7,
8,
9.000000001
]
}
},
{
"name": "svd_rank_one",
"upstream_ref": "source/SMLib/src/SmMatrix.cpp SmMatrix::SVDecompose",
"inputs": {
"op": "svd",
"rows": 3,
"cols": 3,
"m": [
1,
2,
3,
2,
4,
6,
-1,
-2,
-3
]
}
},
{
"name": "svd_huge",
"upstream_ref": "source/SMLib/src/SmMatrix.cpp SmMatrix::SVDecompose",
"inputs": {
"op": "svd",
"rows": 3,
"cols": 3,
"m": [
1e+90,
2e+89,
0,
0,
3e+90,
1e+88,
5e+87,
0,
7e+89
]
}
},
{
"name": "svd_tiny",
"upstream_ref": "source/SMLib/src/SmMatrix.cpp SmMatrix::SVDecompose",
"inputs": {
"op": "svd",
"rows": 3,
"cols": 3,
"m": [
1e-90,
2e-91,
0,
0,
3e-90,
1e-92,
5e-93,
0,
7e-91
]
}
},
{
"name": "svd_tall_4x3",
"upstream_ref": "source/SMLib/src/SmMatrix.cpp SmMatrix::SVDecompose",
"inputs": {
"op": "svd",
"rows": 4,
"cols": 3,
"m": [
1,
2,
3,
4,
5,
6,
7,
8,
9,
10,
11,
12
]
}
},
{
"name": "svd_column_5x1",
"upstream_ref": "source/SMLib/src/SmMatrix.cpp SmMatrix::SVDecompose",
"inputs": {
"op": "svd",
"rows": 5,
"cols": 1,
"m": [
3,
-4,
0,
1,
2
]
}
},
{
"name": "svd_scalar_1x1",
"upstream_ref": "source/SMLib/src/SmMatrix.cpp SmMatrix::SVDecompose",
"inputs": {
"op": "svd",
"rows": 1,
"cols": 1,
"m": [
-2.5
]
}
},
{
"name": "svd_wide_2x3_sm_err",
"upstream_ref": "source/SMLib/src/SmMatrix.cpp SmMatrix::SVDecompose",
"inputs": {
"op": "svd",
"rows": 2,
"cols": 3,
"m": [
1,
2,
3,
4,
5,
6
]
}
}
]
}
File diff suppressed because it is too large Load Diff
+15
View File
@@ -20,6 +20,7 @@ enum ScMatrixOp : int32_t
OP_SOLVE = 8, // SolveLinearSystem(v): out = x
OP_EIGEN = 9, // EigenSystemJacobiSolver: out = values, vectors, matrix; ints = [iter]
OP_EIGEN3 = 10, // Compute3x3EigenVectors: out = values, vectors (rows = vectors), matrix
OP_SVD = 11, // SVDecompose(W, V): out = W, V, matrix after (= U)
};
static void put_matrix(const SmMatrix& m, double* out, int64_t& k)
@@ -209,6 +210,20 @@ SC_EXPORT int32_t sc_geom_matrix_op(int32_t op, int64_t rows, int64_t cols, cons
}
break;
}
case OP_SVD:
{
SmTArray<double> sW;
SmMatrix sVm(1, 1);
*status = sM.SVDecompose(sW, sVm);
if (*status == SM_SUCCESS)
{
for (ULONG r = 0; r < sW.GetSize(); r++)
out[k++] = sW[r];
put_matrix(sVm, out, k);
put_matrix(sM, out, k);
}
break;
}
default:
return SC_ERR_ARG;
}
+15 -5
View File
@@ -67,10 +67,13 @@ void put_geometry(std::vector<double>& o, SmBrep* b) {
// Args: editing, placement O/X/Y (9). Op 0 exports the fixture graph (`0 graph`), op 2 its
// native face and vertex lists (`nf face* nv vertex*`); op 1
// runs Mirror: `record_len status geometry... 0 graph` of the mirrored B-rep.
// runs Mirror: `record_len status geometry... 0 graph` of the mirrored B-rep. Op 3 takes
// 18 args (editing, bScaleTol, the 4x4 matrix row-major) and runs
// SmBrep::Transform(const double[4][4], SmBoolean) with the same output as op 1.
SC_EXPORT int32_t sc_topology_brep_mirror_eval(int32_t op, const char* path, const double* input,
int64_t n, const double* a, int64_t na, double* out, int64_t cap, int64_t* out_len) {
if (op < 0 || op > 2 || !input || n < 1 || !a || na != 10 || !path || !out_len) return SC_ERR_ARG;
if (op < 0 || op > 3 || !input || n < 1 || !a || na != (op == 3 ? 18 : 10) || !path || !out_len)
return SC_ERR_ARG;
SC_TRY
Rd rd{input, 0};
SmStatus status = SM_SUCCESS;
@@ -103,9 +106,16 @@ SC_EXPORT int32_t sc_topology_brep_mirror_eval(int32_t op, const char* path, con
for (ULONG i = 0; i < vertices.GetSize(); ++i) o.push_back((double)g.find(9, vertices[i]));
return sc_copy_out(o.data(), (int64_t)o.size(), out, cap, out_len);
}
SmAxis2Placement plane(SmPoint3d(a[1], a[2], a[3]), SmVector3d(a[4], a[5], a[6]),
SmVector3d(a[7], a[8], a[9]));
status = b->Mirror(plane);
if (op == 3) {
double m[4][4];
for (int r = 0; r < 4; ++r)
for (int c = 0; c < 4; ++c) m[r][c] = a[2 + 4 * r + c];
status = b->Transform(m, a[1] != 0 ? TRUE : FALSE);
} else {
SmAxis2Placement plane(SmPoint3d(a[1], a[2], a[3]), SmVector3d(a[4], a[5], a[6]),
SmVector3d(a[7], a[8], a[9]));
status = b->Mirror(plane);
}
std::vector<double> record = {st(status)};
put_geometry(record, b);
std::vector<double> result = {(double)record.size()};
+1 -1
View File
@@ -9,4 +9,4 @@ tests = ["prog_test:my_test_surface_closest_point", "prog_test:my_test_drop_poin
status = "green"
claimed_by = "swarm-geom.matrix-a1"
claimed_at = "2026-10-12T19:48:03Z"
notes = "20 fns; largest: SmMatrix::EigenSystemJacobiSolver, SmMatrix::LUDecompose, SmMatrix::SolveLinearSystem, SmMatrix::Compute3x3EigenVectors; cycle back-edges (call into later units; stub or co-port): geom.vector3d(2); unit-SCC #22 (220 units); 317 upstream test entry points reach this unit; GREEN: difftest 4044 cases (44 fixture + 4000 generated), 0 deferred, max_rel_err 0 (bit-identical); 100k-seed run clean; back-edges sm_sqrt/SmVector3d::Set (geom.vector3d) co-ported as private helpers"
notes = "20 fns; largest: SmMatrix::EigenSystemJacobiSolver, SmMatrix::LUDecompose, SmMatrix::SolveLinearSystem, SmMatrix::Compute3x3EigenVectors; cycle back-edges (call into later units; stub or co-port): geom.vector3d(2); unit-SCC #22 (220 units); 317 upstream test entry points reach this unit; GREEN: difftest 4044 cases (44 fixture + 4000 generated), 0 deferred, max_rel_err 0 (bit-identical); 100k-seed run clean; back-edges sm_sqrt/SmVector3d::Set (geom.vector3d) co-ported as private helpers SVDecompose and Transpose ported (#1570, PR #1595): Matrix::sv_decompose / Matrix::transpose; difftest 4658 cases (adds 14 svd fixtures + 600 generated svd####), 0 deferred, max_rel_err 0."
+1 -1
View File
@@ -9,4 +9,4 @@ tests = []
status = "green"
claimed_by = ""
claimed_at = ""
notes = "GREEN: difftest 196 cases (68 fixtures: 17 upstream .smb solids x 4 planes, + 128 generated), 196 pass, 0 deferred, 0 failed, max_rel_err=0e0 (--seeds 1000: 1068/1068; --perturb: 196/196 detected; 3 repeat runs identical, ASLR on). NOT IN callgraph/units.json: SmBrep::Mirror and its callees are in harness/callgraph/unused_functions.txt (no upstream test or SM_API entry point reaches them). Owns SmBrep::Mirror (SmBrep.cpp:23449), SmAxis2Placement::MirrorPoint / MirrorVector, and the CreateMirrorCurve / CreateMirrorSurface overrides of SmBSplineSurface, SmOffsetSurface, SmCrvOnSurf and SmCrvInVolume; the SmBSplineCurve and SmSrfInVolume overrides were ported earlier (nurbs.b_spline_curve_p1, surface.srf_in_volume). Rust: topology::brep::mirror::mirror; Axis2Placement::mirror_point / mirror_vector; Curve::create_mirror_curve and Surface::create_mirror_surface dispatch every kind (SM_ERR where upstream's base default applies: Hermite, projected, iso, offset, conic and composite curves; curve bounded and STEP surfaces). The runner compares status, the whole graph and samples of every edge curve, edge use UV trim curve and face surface (type code, domain, points). Uncertified, not reached by a fixture: the OffsetSurface, CrvOnSurf, CrvInVolume and SrfInVolume overrides, the SM_ERR kinds, and bodies with a composite edge (no arena model). Excluded: Solids/fillet_ts.smb (the native vertex list holds a vertex the import does not create) and SmallBox.smb (the runner's native import crashes the process). Deviation: make_missing_nurbs makes the NURBs of analytic surfaces before FindAndSetInfiniteRegion (upstream's bounding box makes them lazily; query.brep--boolean-results-need-their-analytic-nurbs). The runner recomputes Face::shell / Shell::faces on both sides (no upstream counterpart) and sets the native face / vertex list orders on the Rust import."
notes = "GREEN: difftest 196 cases (68 fixtures: 17 upstream .smb solids x 4 planes, + 128 generated), 196 pass, 0 deferred, 0 failed, max_rel_err=0e0 (--seeds 1000: 1068/1068; --perturb: 196/196 detected; 3 repeat runs identical, ASLR on). NOT IN callgraph/units.json: SmBrep::Mirror and its callees are in harness/callgraph/unused_functions.txt (no upstream test or SM_API entry point reaches them). Owns SmBrep::Mirror (SmBrep.cpp:23449), SmAxis2Placement::MirrorPoint / MirrorVector, and the CreateMirrorCurve / CreateMirrorSurface overrides of SmBSplineSurface, SmOffsetSurface, SmCrvOnSurf and SmCrvInVolume; the SmBSplineCurve and SmSrfInVolume overrides were ported earlier (nurbs.b_spline_curve_p1, surface.srf_in_volume). Rust: topology::brep::mirror::mirror; Axis2Placement::mirror_point / mirror_vector; Curve::create_mirror_curve and Surface::create_mirror_surface dispatch every kind (SM_ERR where upstream's base default applies: Hermite, projected, iso, offset, conic and composite curves; curve bounded and STEP surfaces). The runner compares status, the whole graph and samples of every edge curve, edge use UV trim curve and face surface (type code, domain, points). Uncertified, not reached by a fixture: the OffsetSurface, CrvOnSurf, CrvInVolume and SrfInVolume overrides, the SM_ERR kinds, and bodies with a composite edge (no arena model). Excluded: Solids/fillet_ts.smb (the native vertex list holds a vertex the import does not create) and SmallBox.smb (the runner's native import crashes the process). Deviation: make_missing_nurbs makes the NURBs of analytic surfaces before FindAndSetInfiniteRegion (upstream's bounding box makes them lazily; query.brep--boolean-results-need-their-analytic-nurbs). The runner recomputes Face::shell / Shell::faces on both sides (no upstream counterpart) and sets the native face / vertex list orders on the Rust import. SmBrep::Transform(const double[4][4], SmBoolean) ported here (#1570, PR #1595): topology::brep::mirror::transform_matrix; difftest 341 cases (adds 49 transform_matrix fixtures on box/cylinder/sphere/cone + 96 generated mat####), 0 deferred, max_rel_err 0; mat0076 (near-flat sheared revolution) needs the ray-fire GetTree dispatch of PR #1600 and the SmSurfaceEval::EvalNormal fix of PR #1690 (the brute-force UV trim drop's GlobalPointSolve on the flattened face); with both, 341/341 pass."