|
| | ExactBandIntegral (const Eigen::MatrixXd &V, const Eigen::MatrixXi &F) |
| |
| | ExactBandIntegral (const Eigen::MatrixXd &V, const Eigen::MatrixXi &F, const std::vector< int8_t > &face_side, int sigma, size_t value_slots=0, std::shared_ptr< const wmtk::Envelope > envelope=nullptr) |
| |
| double | value (const std::array< Eigen::Vector3d, 4 > &t, double delta, Stats *stats=nullptr) const |
| | J(t) = int_t (d - delta) dV for the cell with corners t (either orientation), restricted to the offset side when constructed with one (class doc); 0 for a flat cell (|6 vol| <= 64 eps L^3, L its longest edge). With stats, adds to it.
|
| |
|
double | value (const std::array< Eigen::Vector3d, 4 > &t, double delta, Partition &keep) const |
| | value(t, delta), t's partition read from or written to keep (Partition).
|
| |
| void | reserve_values (size_t slots) const |
| |
| double | value_grad_hess (const std::array< Eigen::Vector3d, 4 > &t, int kv, double delta, Eigen::Vector3d &grad, Eigen::Matrix3d &hess, Stats *stats=nullptr) const |
| | J and its exact gradient and Hessian with respect to corner kv of t, the other three fixed (exact_mesh.value_grad_hess): Reynolds' transport theorem on each smooth piece, a boundary integral over the three faces of t through x = t[kv], plus, on the Hessian, the terms of piece edges sliding over those faces ("cut edges"; they cancel in pairs between exact pieces and remain on rule-(a) pieces). J here sums the same pieces split into smooth parts, so it equals value() to rounding.
|
| |
|
double | value_grad_hess (const std::array< Eigen::Vector3d, 4 > &t, int kv, double delta, Eigen::Vector3d &grad, Eigen::Matrix3d &hess, Partition &keep) const |
| | value_grad_hess(), t's partition read from or written to keep (Partition).
|
| |
|
void | grad_hess (const std::array< Eigen::Vector3d, 4 > &t, int kv, double delta, Eigen::Vector3d &grad, Eigen::Matrix3d &hess, Partition &keep) const |
| | value_grad_hess()'s gradient and Hessian alone, t's partition read from or written to keep: J's own integration over the elementary pieces is skipped (the smoother reads the value from value() and only the derivatives from here – BandVolumeEnergy3D::eval()). The elementary pieces themselves are read from keep when a value() call with keep captured them for this t (bit for bit) and delta, else built as value_grad_hess() builds them; grad and hess equal value_grad_hess()'s bit for bit either way (exact-band-grad-hess-cache).
|
| |
| double | max_distance (const std::array< Eigen::Vector3d, 3 > &f, MaxInfo *info=nullptr) const |
| | max over the triangle f of d (exact_mesh.max_on_triangle), the settle loop's max test. f is cut by the candidates' local regions as a cell is (the same search, region tests, splits, tolerances and scale rule; f a planar convex ring in place of a polytope), and d is convex wherever one member's distance is d, so on each piece the max is at a vertex of the piece:
|
| |
| double | front_face (const std::array< Eigen::Vector3d, 3 > &f, double delta, FaceParts *keep=nullptr) const |
| | front_face(): int over the triangle f of G.n dA, n = (f1 - f0) x (f2 - f0) / |.| (f ordered with the normal pointing out of the band). front_face_derivs(): for the motion of corner kv (the others fixed), the gradient and the one-vertex Hessian of E's dependence through f: grad = n int_f chi (d - delta) phi dA, hess = n (x) int_f chi phi^2 grad d dA + (int_f chi (d - delta) phi dA / |N|) dN/dx
|
| |
|
void | front_face_derivs (const std::array< Eigen::Vector3d, 3 > &f, double delta, int kv, Eigen::Vector3d &grad, Eigen::Matrix3d &hess, FaceParts *keep=nullptr) const |
| |
|
double | bisector_jump (const std::array< Eigen::Vector3d, 4 > &t, double delta, Partition *keep=nullptr) const |
| | D_t = int_t (d - delta) dV - oint_dt G.n dA of the cell t (either orientation): 0 exactly when one feature labels every part of t on the offset side (one field); otherwise the flux of each labelled part's field through its faces not on dt, outward from the part, a face on one of the label's own region planes skipped (the field is tangent to it). The faces between two parts of one label cancel, so the sum is the jumps of G.n where the label changes: on the bisectors, and on region planes that coincide with one (MEASURED, Python, cube_out tet 2234: the prism plane between two coplanar triangles of a cube face is the bisector of the two perpendicular faces at the cube edge it contains) or nearly so (see the .cpp). t's partition read from or written to keep when given (Partition). Throws on a curved bisector, or as value() if the pieces do not tile t.
|
| |
|
void | prepare_front_flux (double delta) |
| | Builds, from the input alone, where bisectors can be: for every pair of features g, h the polytope L_g cap L_h cap box_g(R) cap box_h(R) (and, with a side, cap both features' offset sides), R = reach = 2 delta, kept where it has volume. A bisector point y of g and h with d(y) <= R lies in it: both are nearest at y, so y is in both regions, on the offset side of both planes, and within d_g(y) = d(y) of a point of g (fact (ii) of the class doc), so in g's box grown by R. Called once, before the integral is shared.
|
| |
|
double | front_flux_reach () const |
| | prepare_front_flux()'s R (0 before it), and the number of polytopes it kept.
|
| |
|
size_t | front_flux_polytopes () const |
| |
|
bool | front_flux_flag (const std::array< Eigen::Vector3d, 4 > &t) const |
| | Whether the cell t can meet a bisector at points with d <= front_flux_reach(): t's box meets a kept polytope's box and no face plane of that polytope has all four corners strictly outside it (beyond rounding). Conservative: a cell flagged without a bisector gets D_t = 0 up to rounding (exactly 0 when one feature labels all of it). Says nothing about points with d > reach: the caller flags a cell whose bound on d (distance_bound()) exceeds the reach itself (BandVolumeEnergy3D).
|
| |
| bool | front_flux_flag (const Eigen::Vector3d *t, int n) const |
| |
|
double | distance (const Eigen::Vector3d &x) const |
| | d(x), the distance to the input: min over the triangles of the point-triangle distance.
|
| |
| double | distance_bound (const std::array< Eigen::Vector3d, 4 > &t) const |
| |
|
double | distance_bound (const Eigen::Vector3d *t, const int n) const |
| | The same bound for the n points t; a subset's bound is never larger.
|
| |
|
size_t | num_faces () const |
| |
|
| void | localize (const Eigen::Vector3d *t, int n, Local &w, Stats *stats) const |
| |
| bool | partition (const std::array< Eigen::Vector3d, 4 > &t, Work &w, Stats *stats) const |
| |
| const Work * | cell_partition (const std::array< Eigen::Vector3d, 4 > &t, Stats *stats, Work *own, Partition *keep) const |
| |
|
double | evaluate (const std::array< Eigen::Vector3d, 4 > &t, double delta, Stats *stats, Partition *keep) const |
| | value() without the value table.
|
| |
|
double | value_tabled (const std::array< Eigen::Vector3d, 4 > &t, double delta, Partition *keep) const |
| | value(t, delta, stats, keep) through the value table (value()).
|
| |
| double | value_grad_hess_impl (const std::array< Eigen::Vector3d, 4 > &t, int kv, double delta, Eigen::Vector3d &grad, Eigen::Matrix3d &hess, Stats *stats, Partition *keep, bool want_value) const |
| |
|
void | tri_partition (const std::array< Eigen::Vector3d, 3 > &f, TriWork &w) const |
| | _tri_partition: the same for the face f, its pieces planar rings.
|
| |
|
double | upper_bound (const Eigen::Vector3d *t, int n) const |
| | U: min over all triangles of max over the n points t of d_P.
|
| |
| void | check_cell_pieces (const std::array< Eigen::Vector3d, 4 > &t, double delta, const Work &w, double vol_sum, double empty_vol) const |
| |
|
void | check_face_pieces (const std::array< Eigen::Vector3d, 3 > &f, const TriWork &w, const char *who) const |
| |
| void | label_face (const std::array< Eigen::Vector3d, 3 > &f, FaceWork &fw) const |
| |
|
const FaceWork & | face_parts (const std::array< Eigen::Vector3d, 3 > &f, FaceParts *keep) const |
| | f's labelled parts, from keep when it holds f's, else made (into keep, or a scratch).
|
| |
|
int | build_ff_node (int begin, int end) |
| |
|
void | candidates (const Eigen::Vector3d &lo, const Eigen::Vector3d &hi, double r, std::vector< int > &out) const |
| | Feature indices whose box is within r of [lo, hi], ascending.
|
| |
The band-volume integral of one cell, in closed form:
J(t) = int_t (d(x) - delta) dV, d = Euclidean distance to the input triangle mesh.
A port of the VALUE path of the Python reference exact_mesh.integrate_mesh (and of the exact_tri closed forms it calls), scratchpad exact_integral/ of 2026-10-10/11: same algorithm, same switch constants, same operation order, so it reproduces the reference's numbers to rounding (tests/test_exact_band_integral.cpp). Used by EXPERIMENTAL_band_volume_rule "exact" (the band cell term) and, through max_distance(), by the settle loop's max test.
Features and local regions. The input's features are its face interiors, edges and vertices (duplicate vertices merged). The local region L_g of a feature g is an intersection of half-spaces built from its incident features only: a face's prism (inside the planes through its edges perpendicular to it), an edge's slab between the planes through its ends perpendicular to it, cut by the plane through the edge perpendicular to each incident face (away from that face's third vertex), a vertex's cone {x : (x - v).(w - v) <= 0 for every edge (v, w)}. Two facts make d = min over {g : x in L_g} of d_g exact, d_g the plane, line or point distance of g: (i) the feature holding x's nearest point in its relative interior has x in its region, and d_g(x) = d(x); (ii) on L_g, d_g is the distance to an actual point of g, so d_g >= d.
Per cell: 0. Candidates: U = min over triangles P of max over t's corners of d_P is an upper bound of d on t (d_P convex), so every feature whose bounding box is within U of t's box is kept. U is the exact min over ALL triangles (a box search); the Python takes the min over the 32 triangles nearest t by box distance, MEASURED equal on all 14180 cells checked (the five stored cases, sample + 3000 random each), although ties at rank 32 occur in up to 634 of 3240 cells, where numpy's sort order decides which 32 it takes.
- t is cut by the candidates' local regions into convex pieces of constant member set S = {g : piece inside L_g}. Members that cannot attain the minimum on a piece are dropped (lower bound >= another member's upper bound); this changes no value.
- |S| = 1: the closed form of d_g. S all faces: split by the planes and the bisectors of the (then affine) distances, exact. S all vertices: split by the bisector planes, exact.
- Otherwise ("rule (a)") the boundary between members is curved and the piece returns min over g in S of int_piece d_g >= int_piece d: an UPPER BOUND, not the integral. Stats::rule_a_vol reports its volume (0.15-0.35% of the band volume on t100026).
One-feature integrals (exact_tri.py): the face's |n.(x - a)| is affine on each side of its plane (volume times centroid value); the line and point distances by Euler's theorem for degree-1 homogeneous integrands, int_Q f = 1/4 sum_F h_F int_F f dA, faces by the divergence theorem in their plane, edge integrals in sqrt, asinh and atan. Two pieces are convergent series with explicit truncation bounds, not closed forms: a face nearly parallel to an edge line (eps_F < 0.1) uses a Taylor series in the variation of the line's height over the face, per edge either along the edge (ratio <= 1/4, at most 29 terms) or in eta (ratio < 1/20, at most 13 terms), each stopped at its edge's actual ratio with the same tail bound.
Conditioning (Python, test_exact_tri.py "far"): a piece of size L at distance D from the singular point or line loses about eps (D/L)^2 in the Euler sums.
The offset side (constructed with a side): the band B is the region between the input and the front, so a point on the far side of the input surface from the offset is not in it, and value() returns J(t) = int_{t cap B} (d - delta) dV. The side of x is decided by its nearest feature g, by the sign of (x - c).N_g for its nearest point c (the pseudonormal test, Baerentzen-Aanaes 2005): N_g = the face normal (face), n1 + n2 (edge), the angle-weighted pseudonormal (vertex). On g's local region c is g's own nearest point to x, so the test is the fixed plane {(x - q_g).N_g = 0}, q_g a point of g:
- a one-member piece is clipped by its member's plane and only the offset-side part counts; a faces-only sub-piece (one face nearest, on one side of its plane) is kept or dropped whole; a vertices-only sub-piece is clipped by its vertex's plane;
- a rule-(a) piece is split by its members' planes. Where all planes give the same side, the true side is that one (it is the nearest member's): the part is dropped, or takes rule (a)'s min. Where they disagree (D) the true side is still one for all of D: the plane test of a feature nearest at x gives x's true side, so members of opposite verdicts are never nearest together in D. It is the side of the member nearest at D's centroid: D is dropped, or takes rule (a)'s min over the members that put it on the offset side, the others being nearest nowhere in it. Until 2026-10-14 D took vol(D) max(0, min over members of (max over D's vertices of d_g) - delta), an upper bound that depended on how the cell cut D: MEASURED on 100033 (plates 2.9 thick, delta 2.49), up to 1.24 between two triangulations of one region, and a swap pass that cycled. The plain min over ALL members' own clipped integrals is no bound: MEASURED, two plates 0.02 apart, delta 0.2, a cell wholly inside the lower plate (true value 0): the upper plate's rim edge, whose region held the cell, put it on the offset side, and that min was -1.32e-5 (the lower plate's face, nearest there, drops it).
- -delta vol(t) becomes -delta times the kept volume; a cell where nothing is clipped keeps the unclipped sum, int d - delta vol(t), bit for bit. Features of isolated triangles (no side), edges without exactly two faces and zero pseudonormals have no plane: every point counts there. Why (exact_mesh.band_integral): the input surface inside the tet mesh may move within its envelope, and just inside the input body the integrand is about -1, so pushing the band's inner boundary into the body lowered E; operations accepted when they lower E ratcheted it inward. MEASURED: on the outward sphere 446 of 881 inner-boundary vertices lay inside the input by turn 12, E_band ended 58.6 bars below its exact minimum (Steiner's formula), and the settle test could not pass.
Thread safety: value() and value_grad_hess() are const. Their state is a per-thread scratch partition and value()'s table of results, whose slots are read and written under a lock (value()); concurrent calls are safe. A Partition (below) is the caller's, used by one thread.
| double wmtk::components::topological_offset::ExactBandIntegral::value |
( |
const std::array< Eigen::Vector3d, 4 > & |
t, |
|
|
double |
delta, |
|
|
Stats * |
stats = nullptr |
|
) |
| const |
J(t) = int_t (d - delta) dV for the cell with corners t (either orientation), restricted to the offset side when constructed with one (class doc); 0 for a flat cell (|6 vol| <= 64 eps L^3, L its longest edge). With stats, adds to it.
The value table (constructed with value_slots > 0; not used with stats, whose counters count every evaluation): J is a function of its 13 inputs alone – t's 12 coordinates in their order, delta – and of the input triangles, fixed. A slot, chosen by a hash of the inputs, holds the last result with its inputs; a call whose inputs equal a slot's, bit for bit, returns its result, which is what the evaluation would return again. Any other call evaluates and takes the slot over: a cell that moved or two cells sharing a slot cost an evaluation, never a different value.
Why: the energy rule reads the band cells of an operation's ring before it, and an unchanged cell is read again by every later operation whose ring holds it. MEASURED (Geneva base, inward, exact rule, the turn-1 collapse group, every evaluation of value() logged by its inputs): the collapse before-sums made 8.5M calls (6.2 band cells per sum), every one of them a repeat of an earlier call, 83% of the group's value() time; the after-sums 0.69M, 62% repeats (a collapse refused, retried in a later round). A table taking every result, of 2^18 / 2^20 / 2^22 slots, would have answered 98.2 / 99.5 / 99.9% of the before-sums and 52 / 59 / 61% of the after-sums.
