From e1e907dff445eb5adeff95e1087b17ae7f906b7e Mon Sep 17 00:00:00 2001 From: xywei Date: Fri, 12 Jun 2026 17:58:04 -0500 Subject: [PATCH 01/15] feat(evals): add near-field template experiment Add an analytic folded fan near-field experiment for issue 61. The prototype derives and documents the mapped point-target and bilinear self-interaction forms, frames the reusable correction hypothesis as a functional expansion in the smooth metric field M(z), and adds a script wrapper for running the baseline check. Tests cover fan-map invariants, log-kernel scale behavior, point-target convergence, and the diagonal metric singular split.\n\nValidation: uv run --extra dev python -m pytest tests/evals/test_nearfield_templates.py tests/test_script_wrappers.py; uv run --extra dev python scripts/run_nearfield_template_experiment.py --order 6; make dev Entire-Checkpoint: 3d4ca5797c60 --- ...2-boxcode-nearfield-template-experiment.md | 71 ++++ docs/index.md | 1 + docs/nearfield-template-experiments.md | 163 ++++++++ scripts/run_nearfield_template_experiment.py | 63 +++ src/cutkit/evals/__init__.py | 22 ++ src/cutkit/evals/nearfield_templates.py | 360 ++++++++++++++++++ tests/evals/test_nearfield_templates.py | 105 +++++ 7 files changed, 785 insertions(+) create mode 100644 docs/exec-plans/completed/20260612-boxcode-nearfield-template-experiment.md create mode 100644 docs/nearfield-template-experiments.md create mode 100755 scripts/run_nearfield_template_experiment.py create mode 100644 src/cutkit/evals/nearfield_templates.py create mode 100644 tests/evals/test_nearfield_templates.py diff --git a/docs/exec-plans/completed/20260612-boxcode-nearfield-template-experiment.md b/docs/exec-plans/completed/20260612-boxcode-nearfield-template-experiment.md new file mode 100644 index 0000000..5c85bee --- /dev/null +++ b/docs/exec-plans/completed/20260612-boxcode-nearfield-template-experiment.md @@ -0,0 +1,71 @@ +# Box-Code Near-Field Template Experiment + +## Objective + +Evaluate whether folded-decomposition charts can support reusable near-field +template corrections for box-code and Volumential-style workflows. + +## Scope + +- Derive the mapped-kernel form for a 2D fan folded piece. +- Include both point-target and target-piece bilinear near-field forms. +- Identify singular factors that can be separated from smooth geometry payloads. +- Frame the reuse hypothesis as a functional expansion in smooth metric fields + `M(z)`, not just as raw runtime geometry data. +- Add a small analytic experiment driver with deterministic tests. +- Document feasibility limits and the next CUTKIT/Volumential boundary. + +## Non-Goals + +- No production Volumential integration in this branch. +- No CAD dependency for the first prototype. +- No claim that arbitrary cut shapes have finite exact correction tables. + +## Acceptance Criteria + +- The derivation states the mapped self-interaction form and its split. +- The derivation records the point-target form, bilinear form, interaction + taxonomy, and measurement criteria from issue 61. +- The prototype records a self or adjacent template experiment with a direct + high-order reference value. +- Tests check the experiment against an analytic or independently computed + invariant. +- Documentation explains whether the idea is feasible and what data must remain + runtime geometry payload. + +## Completed Checklist + +- [x] Read issue 61 and repository orientation docs. +- [x] Work through the 2D fan-map singularity derivation. +- [x] Add an experiment module or script for the near-field template prototype. +- [x] Add tests for the prototype invariants. +- [x] Update docs and doc index if new documentation is added. +- [x] Run targeted tests and `make dev`. +- [x] Move this plan to completed with outcomes. + +## Outcomes + +- The near-field template idea is feasible as a CUTKIT experiment: singular and + nearly singular folded-piece interactions can be expressed on fixed template + domains. +- The reusable singular part is a template-space log singular model. The metric + field `M(z) = DT(z)^T DT(z)` and higher-order terms are smooth geometry data. +- The practical reuse hypothesis is that a functional expansion in `M(z)` and + smooth-remainder data covers the folded-decomposition cases with few modes. +- The first prototype covers analytic straight fan maps, point-target + integrals, bilinear self interactions, a log-kernel scale-law check, and a + diagonal metric-remainder sample. + +## Validation + +- `uv run --extra dev python -m pytest tests/evals/test_nearfield_templates.py tests/test_script_wrappers.py` +- `uv run --extra dev python scripts/run_nearfield_template_experiment.py --order 6` +- `make dev` + +## Remaining Risks + +- Template corrections may be reusable only after parameterized compression, not + as shape-independent lookup tables. +- Self-interaction singular quadrature may need specialized rules beyond the + first direct high-order reference prototype. +- Curved analytic pieces and adjacent-piece cases still need follow-up coverage. diff --git a/docs/index.md b/docs/index.md index ff80f0e..b6fe266 100644 --- a/docs/index.md +++ b/docs/index.md @@ -23,6 +23,7 @@ Repository-local documentation is the system of record for CUTKIT. - `docs/cad-box-batch-interface.md`: consumer-facing CAD object-or-arrays clip/integrate interface - `docs/volumential-handoff.md`: CUTKIT far/near primitive handoff contract for volumential workflows - `docs/meshmode-cut-overlay.md`: meshmode cut-overlay contract and validation semantics +- `docs/nearfield-template-experiments.md`: folded fan near-field template derivation and prototype - `docs/poisson-benchmarks.md`: Poisson-oriented benchmark usage and profiles - `docs/poisson-galerkin-solver.md`: immersed Poisson Galerkin solve + validation workflow - `docs/formdsl-parity-benchmarks.md`: shared IGA + DG-SEM formdsl parity benchmark usage diff --git a/docs/nearfield-template-experiments.md b/docs/nearfield-template-experiments.md new file mode 100644 index 0000000..070db97 --- /dev/null +++ b/docs/nearfield-template-experiments.md @@ -0,0 +1,163 @@ +# Near-Field Template Experiments + +This note records the first mathematical check for issue 61: whether folded +decomposition charts can support a reusable near-field template story for later +box-code or Volumential work. + +## Fan Map + +For one straight-edge folded 2D piece, use the fan chart + +```text +T(r, t) = (1-r) V + r C(t), +C(t) = (1-t) P0 + t P1, +0 <= r,t <= 1. +``` + +Its Jacobian is + +```text +J(r, t) = r det(C(t)-V, C'(t)) = r det(P0-V, P1-V). +``` + +The factor `r` is the Duffy-style apex degeneracy already used by +`cutkit.quadrature.folded2d`. It is not a new singularity for physical +integration; it cancels area at the apex. + +## Mapped Kernel + +For the 2D Laplace kernel + +```text +G(x, y) = -log(|x-y|)/(2*pi), +``` + +the point-target near-field integral is + +```text +u(x) = integral_[0,1]^2 G(x, T_e(r,t)) rho(T_e(r,t)) J_e(r,t) dr dt. +``` + +This is the cheapest diagnostic path because it isolates source-piece mapping, +source density, and target location before introducing target basis functions. +For target pieces with a second map `T_f(eta)`, the bilinear template integral is + +```text +A_fe = integral_[0,1]^2 integral_[0,1]^2 + psi_f(eta) G(T_f(eta), T_e(xi)) phi_e(xi) + J_f(eta) J_e(xi) dxi deta. +``` + +For the self-piece case, this specializes to + +```text +I = integral_[0,1]^2 integral_[0,1]^2 + G(T(r,t), T(r',t')) J(r,t) J(r',t') dr dt dr' dt'. +``` + +The experiment should classify interactions as: + +- self piece, `e = f`; +- edge-adjacent pieces sharing a boundary segment; +- vertex or seed-adjacent pieces; +- near but disjoint pieces; +- well-separated pieces, where ordinary tensor-product quadrature should already work. + +## Singular Split + +Near a non-apex diagonal point `z=(r,t)`, write `delta = z' - z`. Then + +```text +T(z + delta) - T(z) = DT(z) delta + O(|delta|^2), +|T(z + delta) - T(z)| = sqrt(delta^T M(z) delta) (1 + O(|delta|)), +M(z) = DT(z)^T DT(z). +``` + +Therefore + +```text +G(T(z), T(z+delta)) += -log(sqrt(delta^T M(z) delta))/(2*pi) + smooth remainder. +``` + +The reusable part would be template-space singular model integrals or correction +operators. The runtime payload remains map coefficients, Jacobian data, source +coefficients, target metadata, and smooth-remainder moment/interpolation data. + +## Feasibility Result + +The answer is qualified yes. Singular or nearly singular folded-piece +interactions can be moved to fixed template domains, but not to a finite exact +table independent of geometry. The singular coordinate form is reusable; the +metric and higher-order geometry terms enter through smooth geometry-dependent +payloads. + +The practical hypothesis is stronger than carrying `M(z)` pointwise at runtime: +because the folded maps have smooth, low-parameter structure away from collapsed +seed faces, the metric field `M(z)` should vary smoothly over the template and +across the geometry families produced by folded decomposition. A functional +expansion in the metric data, for example moments, interpolation coefficients, +or a low-rank basis in `M` and the smooth remainder, may cover enough practical +cut-piece cases to make reusable near-field corrections worthwhile. + +The open numerical question is therefore whether the smooth dependence on `M`, +seed location, and curve coefficients is compact enough to approximate by such a +functional expansion for useful geometry families. + +## Measurements + +For each geometry and interaction case, record: + +- reference value; +- ordinary pulled-forward tensor-product quadrature error; +- template singular-correction error; +- smooth-remainder approximation error; +- dependence on quadrature order; +- dependence on seed location and curve coefficients; +- size, rank, and smoothness of geometry-parameterized correction data; +- how many metric-field expansion modes are needed for the observed folded + decomposition cases. + +The first geometry family is CAD-independent: straight fan pieces, then +quadratic Bezier arcs and mildly curved cubics, with seed locations chosen to +produce positive, negative, and folded orientations. OpenCascade should only be +used later as a stress-test backend. + +## Prototype + +`cutkit.evals.nearfield_templates` implements: + +- `FanTemplateMap2D` for the analytic fan map. +- `point_target_laplace_potential(...)` for point-target source integrals with + polynomial template densities. +- `self_interaction_laplace(...)` for a direct high-order mapped self integral. +- `diagonal_remainder_sample(...)` for the metric singular split. +- `run_nearfield_template_experiment(...)` for the baseline feasibility check. + +The script wrapper is: + +```bash +uv run python scripts/run_nearfield_template_experiment.py --order 12 +``` + +The experiment checks the exact scale law for the 2D log kernel. If the fan is +scaled by `lambda`, then + +```text +I_lambda = lambda^4 (I - log(lambda) A^2/(2*pi)), +``` + +where `A` is the signed area of the unscaled fan piece. It also records diagonal +remainder samples that shrink as the source/target template coordinates coalesce +and a point-target low-order-vs-reference error. + +## Next Decision + +The idea is feasible as a CUTKIT experiment, with these limits: + +- Far-field source clouds can remain ordinary signed quadrature sources. +- Near-field correction reuse should target template singular bases plus + functional expansions in the smooth metric field `M(z)` and higher geometry + terms. +- Volumential should still own tree/list composition; CUTKIT should export the + local geometry/operator payloads needed by those lists. diff --git a/scripts/run_nearfield_template_experiment.py b/scripts/run_nearfield_template_experiment.py new file mode 100755 index 0000000..1e4442f --- /dev/null +++ b/scripts/run_nearfield_template_experiment.py @@ -0,0 +1,63 @@ +#!/usr/bin/env python3 +"""Run the folded fan near-field template experiment.""" + +from __future__ import annotations + +import argparse + +from cutkit.evals import run_nearfield_template_experiment + + +def _parse_args() -> argparse.Namespace: + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument( + "--order", + type=int, + default=12, + help="Gauss-Legendre order for the target template coordinates.", + ) + parser.add_argument( + "--scale-factor", + type=float, + default=1.75, + help="Positive geometry scale factor for the log-kernel scale-law check.", + ) + return parser.parse_args() + + +def main() -> int: + args = _parse_args() + result = run_nearfield_template_experiment( + order=args.order, + scale_factor=args.scale_factor, + ) + + print("field | value") + print("--- | ---") + print(f"order | {result.order}") + print(f"source_order | {result.source_order}") + print(f"near_point_target | {result.near_point_target}") + print(f"point_target_reference | {result.point_target_reference:.16e}") + print(f"point_target_low_order | {result.point_target_low_order:.16e}") + print(f"point_target_abs_error | {result.point_target_abs_error:.16e}") + print(f"signed_area | {result.signed_area:.16e}") + print(f"self_interaction | {result.self_interaction:.16e}") + print(f"scale_factor | {result.scale_factor:.16e}") + print(f"scaled_self_interaction | {result.scaled_self_interaction:.16e}") + print( + f"expected_scaled_self_interaction | {result.expected_scaled_self_interaction:.16e}" + ) + print(f"scaled_abs_error | {result.scaled_abs_error:.16e}") + print("") + print("delta | physical_distance | model_distance | remainder") + print("--- | --- | --- | ---") + for sample in result.diagonal_remainders: + print( + f"{sample.delta:.1e} | {sample.physical_distance:.16e} | " + f"{sample.model_distance:.16e} | {sample.remainder:.16e}" + ) + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/src/cutkit/evals/__init__.py b/src/cutkit/evals/__init__.py index 110c262..245c6d2 100644 --- a/src/cutkit/evals/__init__.py +++ b/src/cutkit/evals/__init__.py @@ -74,6 +74,18 @@ run_poisson_galerkin_benchmark, solve_trimmed_poisson_galerkin, ) +from cutkit.evals.nearfield_templates import ( + DiagonalRemainderSample, + FanTemplateMap2D, + NearfieldTemplateExperiment, + diagonal_remainder_sample, + expected_scaled_laplace_self_interaction, + laplace_log_kernel, + metric_model_distance, + point_target_laplace_potential, + run_nearfield_template_experiment, + self_interaction_laplace, +) __all__ = [ "case_to_fixture_payload", @@ -138,6 +150,16 @@ "build_poisson_galerkin_geometry_snapshot", "run_poisson_galerkin_benchmark", "solve_trimmed_poisson_galerkin", + "DiagonalRemainderSample", + "FanTemplateMap2D", + "NearfieldTemplateExperiment", + "diagonal_remainder_sample", + "expected_scaled_laplace_self_interaction", + "laplace_log_kernel", + "metric_model_distance", + "point_target_laplace_potential", + "run_nearfield_template_experiment", + "self_interaction_laplace", "FormDslParityBenchmark", "FormDslParityRow", "FormDslMultipatchStressBenchmark", diff --git a/src/cutkit/evals/nearfield_templates.py b/src/cutkit/evals/nearfield_templates.py new file mode 100644 index 0000000..d28e5d7 --- /dev/null +++ b/src/cutkit/evals/nearfield_templates.py @@ -0,0 +1,360 @@ +"""Near-field template experiments for folded 2D fan pieces.""" + +from __future__ import annotations + +from dataclasses import dataclass +from math import hypot, log, pi, sqrt +from typing import Callable + +from cutkit.quadrature import gauss_legendre_01 + +Point2D = tuple[float, float] +TemplateDensity2D = Callable[[float, float], float] + + +def unit_template_density(_r: float, _t: float) -> float: + """Return the constant template density used by default experiments.""" + + return 1.0 + + +def _sub(lhs: Point2D, rhs: Point2D) -> Point2D: + return (lhs[0] - rhs[0], lhs[1] - rhs[1]) + + +def _dot(lhs: Point2D, rhs: Point2D) -> float: + return lhs[0] * rhs[0] + lhs[1] * rhs[1] + + +def _cross(lhs: Point2D, rhs: Point2D) -> float: + return lhs[0] * rhs[1] - lhs[1] * rhs[0] + + +def _scale_point(point: Point2D, scale: float) -> Point2D: + return (scale * point[0], scale * point[1]) + + +@dataclass(frozen=True) +class FanTemplateMap2D: + """Straight-edge folded fan chart ``T(r, t) = (1-r)V + r C(t)``.""" + + vertex: Point2D + edge_start: Point2D + edge_end: Point2D + + def point(self, r: float, t: float) -> Point2D: + """Map ``(r, t)`` in ``[0, 1]^2`` to physical space.""" + + one_minus_t = 1.0 - t + curve_x = one_minus_t * self.edge_start[0] + t * self.edge_end[0] + curve_y = one_minus_t * self.edge_start[1] + t * self.edge_end[1] + return ( + (1.0 - r) * self.vertex[0] + r * curve_x, + (1.0 - r) * self.vertex[1] + r * curve_y, + ) + + @property + def boundary_det(self) -> float: + """Return ``det(C(t)-V, C'(t))``, constant for a straight edge.""" + + return _cross( + _sub(self.edge_start, self.vertex), + _sub(self.edge_end, self.edge_start), + ) + + @property + def signed_area(self) -> float: + """Return the oriented physical area of the fan piece.""" + + return 0.5 * self.boundary_det + + def signed_jacobian(self, r: float) -> float: + """Return the signed Jacobian of the ``(r, t)`` fan chart.""" + + return r * self.boundary_det + + def local_metric( + self, + r: float, + t: float, + ) -> tuple[tuple[float, float], tuple[float, float]]: + """Return ``DT(r,t)^T DT(r,t)`` for the fan chart.""" + + radial = _sub(self.point(1.0, t), self.vertex) + tangent_base = _sub(self.edge_end, self.edge_start) + tangent = (r * tangent_base[0], r * tangent_base[1]) + return ( + (_dot(radial, radial), _dot(radial, tangent)), + (_dot(tangent, radial), _dot(tangent, tangent)), + ) + + def scaled(self, factor: float) -> FanTemplateMap2D: + """Return the same template scaled about the origin.""" + + return FanTemplateMap2D( + vertex=_scale_point(self.vertex, factor), + edge_start=_scale_point(self.edge_start, factor), + edge_end=_scale_point(self.edge_end, factor), + ) + + +@dataclass(frozen=True) +class DiagonalRemainderSample: + """One mapped-kernel singular split sample near the diagonal.""" + + delta: float + physical_distance: float + model_distance: float + laplace_kernel: float + metric_model_kernel: float + remainder: float + + +@dataclass(frozen=True) +class NearfieldTemplateExperiment: + """Summary of one fan-chart self-interaction experiment.""" + + fan: FanTemplateMap2D + order: int + source_order: int + near_point_target: Point2D + point_target_reference: float + point_target_low_order: float + point_target_abs_error: float + self_interaction: float + signed_area: float + scale_factor: float + scaled_self_interaction: float + expected_scaled_self_interaction: float + scaled_abs_error: float + diagonal_remainders: tuple[DiagonalRemainderSample, ...] + + +def laplace_log_kernel(target: Point2D, source: Point2D) -> float: + """Return the 2D Laplace fundamental solution ``-log(|x-y|)/(2*pi)``.""" + + distance = hypot(target[0] - source[0], target[1] - source[1]) + if distance <= 0.0: + raise ValueError("Laplace log kernel is singular at coincident points") + return -log(distance) / (2.0 * pi) + + +def metric_model_distance( + fan: FanTemplateMap2D, + *, + r: float, + t: float, + delta_r: float, + delta_t: float, +) -> float: + """Return the local metric distance for a template-space displacement.""" + + metric = fan.local_metric(r, t) + distance_squared = ( + metric[0][0] * delta_r * delta_r + + 2.0 * metric[0][1] * delta_r * delta_t + + metric[1][1] * delta_t * delta_t + ) + if distance_squared <= 0.0: + raise ValueError("metric model distance must be positive") + return sqrt(distance_squared) + + +def self_interaction_laplace( + fan: FanTemplateMap2D, + *, + order: int, + source_order: int | None = None, + source_density: TemplateDensity2D | None = None, + target_density: TemplateDensity2D | None = None, +) -> float: + """Approximate the mapped self interaction on one fan chart.""" + + if order < 1: + raise ValueError("order must be positive") + if source_order is None: + source_order = order + 1 + if source_order < 1: + raise ValueError("source_order must be positive") + if source_density is None: + source_density = unit_template_density + if target_density is None: + target_density = unit_template_density + + target_nodes, target_weights = gauss_legendre_01(order) + source_nodes, source_weights = gauss_legendre_01(source_order) + total = 0.0 + + for target_r, target_wr in zip(target_nodes, target_weights, strict=True): + target_j = fan.signed_jacobian(target_r) + for target_t, target_wt in zip(target_nodes, target_weights, strict=True): + target = fan.point(target_r, target_t) + target_weight = ( + target_density(target_r, target_t) * target_j * target_wr * target_wt + ) + for source_r, source_wr in zip(source_nodes, source_weights, strict=True): + source_j = fan.signed_jacobian(source_r) + for source_t, source_wt in zip( + source_nodes, source_weights, strict=True + ): + source = fan.point(source_r, source_t) + source_weight = ( + source_density(source_r, source_t) + * source_j + * source_wr + * source_wt + ) + total += ( + laplace_log_kernel(target, source) + * target_weight + * source_weight + ) + + return total + + +def point_target_laplace_potential( + fan: FanTemplateMap2D, + target: Point2D, + *, + order: int, + source_density: TemplateDensity2D | None = None, +) -> float: + """Approximate a folded-piece source integral at one point target.""" + + if order < 1: + raise ValueError("order must be positive") + if source_density is None: + source_density = unit_template_density + + nodes, weights = gauss_legendre_01(order) + total = 0.0 + for r, wr in zip(nodes, weights, strict=True): + jacobian = fan.signed_jacobian(r) + for t, wt in zip(nodes, weights, strict=True): + source = fan.point(r, t) + total += ( + laplace_log_kernel(target, source) + * source_density(r, t) + * jacobian + * wr + * wt + ) + return total + + +def expected_scaled_laplace_self_interaction( + base_interaction: float, + signed_area: float, + scale_factor: float, +) -> float: + """Return the exact 2D log-kernel scale law for a scaled fan piece.""" + + if scale_factor <= 0.0: + raise ValueError("scale_factor must be positive") + return scale_factor**4 * ( + base_interaction - log(scale_factor) * signed_area * signed_area / (2.0 * pi) + ) + + +def diagonal_remainder_sample( + fan: FanTemplateMap2D, + *, + r: float, + t: float, + delta: float, + direction: tuple[float, float] = (1.0, 0.5), +) -> DiagonalRemainderSample: + """Compare the full mapped kernel with the local metric singular model.""" + + if delta <= 0.0: + raise ValueError("delta must be positive") + delta_r = delta * direction[0] + delta_t = delta * direction[1] + source = fan.point(r, t) + target = fan.point(r + delta_r, t + delta_t) + physical_distance = hypot(target[0] - source[0], target[1] - source[1]) + model_distance = metric_model_distance( + fan, + r=r, + t=t, + delta_r=delta_r, + delta_t=delta_t, + ) + laplace_kernel = -log(physical_distance) / (2.0 * pi) + metric_model_kernel = -log(model_distance) / (2.0 * pi) + return DiagonalRemainderSample( + delta=delta, + physical_distance=physical_distance, + model_distance=model_distance, + laplace_kernel=laplace_kernel, + metric_model_kernel=metric_model_kernel, + remainder=laplace_kernel - metric_model_kernel, + ) + + +def run_nearfield_template_experiment( + *, + order: int = 12, + scale_factor: float = 1.75, +) -> NearfieldTemplateExperiment: + """Run the baseline analytic fan near-field template experiment.""" + + fan = FanTemplateMap2D( + vertex=(0.0, 0.0), + edge_start=(1.0, 0.0), + edge_end=(0.35, 0.9), + ) + source_order = order + 1 + self_value = self_interaction_laplace( + fan, + order=order, + source_order=source_order, + ) + near_point_target = (0.47, 0.42) + + def point_density(r: float, t: float) -> float: + return 1.0 + r - 0.5 * t + + point_reference = point_target_laplace_potential( + fan, + near_point_target, + order=max(order + 8, 16), + source_density=point_density, + ) + point_low_order = point_target_laplace_potential( + fan, + near_point_target, + order=max(order // 2, 2), + source_density=point_density, + ) + scaled_value = self_interaction_laplace( + fan.scaled(scale_factor), + order=order, + source_order=source_order, + ) + expected_scaled = expected_scaled_laplace_self_interaction( + self_value, + fan.signed_area, + scale_factor, + ) + remainders = tuple( + diagonal_remainder_sample(fan, r=0.43, t=0.37, delta=delta) + for delta in (1.0e-1, 1.0e-2, 1.0e-3) + ) + return NearfieldTemplateExperiment( + fan=fan, + order=order, + source_order=source_order, + near_point_target=near_point_target, + point_target_reference=point_reference, + point_target_low_order=point_low_order, + point_target_abs_error=abs(point_low_order - point_reference), + self_interaction=self_value, + signed_area=fan.signed_area, + scale_factor=scale_factor, + scaled_self_interaction=scaled_value, + expected_scaled_self_interaction=expected_scaled, + scaled_abs_error=abs(scaled_value - expected_scaled), + diagonal_remainders=remainders, + ) diff --git a/tests/evals/test_nearfield_templates.py b/tests/evals/test_nearfield_templates.py new file mode 100644 index 0000000..23cb351 --- /dev/null +++ b/tests/evals/test_nearfield_templates.py @@ -0,0 +1,105 @@ +from __future__ import annotations + +import pytest + +from cutkit.evals import ( + FanTemplateMap2D, + diagonal_remainder_sample, + expected_scaled_laplace_self_interaction, + point_target_laplace_potential, + run_nearfield_template_experiment, + self_interaction_laplace, +) + + +def test_fan_map_area_and_metric_are_consistent() -> None: + fan = FanTemplateMap2D( + vertex=(0.0, 0.0), + edge_start=(1.0, 0.0), + edge_end=(0.25, 0.5), + ) + + assert fan.signed_area == pytest.approx(0.25) + metric = fan.local_metric(0.5, 0.25) + assert metric[0][0] > 0.0 + assert metric[1][1] > 0.0 + assert metric[0][1] == pytest.approx(metric[1][0]) + + +def test_self_interaction_obeys_log_kernel_scale_law() -> None: + fan = FanTemplateMap2D( + vertex=(0.0, 0.0), + edge_start=(1.0, 0.0), + edge_end=(0.35, 0.9), + ) + scale_factor = 2.25 + value = self_interaction_laplace(fan, order=5, source_order=6) + scaled = self_interaction_laplace( + fan.scaled(scale_factor), + order=5, + source_order=6, + ) + expected = expected_scaled_laplace_self_interaction( + value, + fan.signed_area, + scale_factor, + ) + + assert scaled == pytest.approx(expected, abs=1.0e-13) + + +def test_metric_singular_remainder_shrinks_toward_diagonal() -> None: + fan = FanTemplateMap2D( + vertex=(0.0, 0.0), + edge_start=(1.0, 0.0), + edge_end=(0.35, 0.9), + ) + coarse = diagonal_remainder_sample(fan, r=0.43, t=0.37, delta=1.0e-2) + fine = diagonal_remainder_sample(fan, r=0.43, t=0.37, delta=1.0e-4) + + assert abs(fine.remainder) < abs(coarse.remainder) + assert fine.physical_distance == pytest.approx(fine.model_distance, rel=1.0e-3) + + +def test_point_target_path_converges_for_near_disjoint_target() -> None: + fan = FanTemplateMap2D( + vertex=(0.0, 0.0), + edge_start=(1.0, 0.0), + edge_end=(0.35, 0.9), + ) + target = (0.8, 0.7) + + def density(r: float, t: float) -> float: + return 1.0 + r * t + + coarse = point_target_laplace_potential( + fan, + target, + order=4, + source_density=density, + ) + fine = point_target_laplace_potential( + fan, + target, + order=10, + source_density=density, + ) + reference = point_target_laplace_potential( + fan, + target, + order=18, + source_density=density, + ) + + assert abs(fine - reference) < abs(coarse - reference) + + +def test_baseline_experiment_records_feasibility_checks() -> None: + result = run_nearfield_template_experiment(order=5, scale_factor=1.5) + + assert result.scaled_abs_error < 1.0e-13 + assert result.point_target_abs_error > 0.0 + assert result.diagonal_remainders[-1].delta == pytest.approx(1.0e-3) + assert abs(result.diagonal_remainders[-1].remainder) < abs( + result.diagonal_remainders[0].remainder + ) From e987ed6079ae5482d1a9d7e2f7ee3e4fb09efdf8 Mon Sep 17 00:00:00 2001 From: xywei Date: Fri, 12 Jun 2026 18:13:17 -0500 Subject: [PATCH 02/15] docs: derive metric expansion tables Document how smooth metric fields in folded chart singular integrals can be expanded around a reference metric to recover fixed template tables times per-chart runtime coefficients. The derivation ties the precomputed-table technique to the fan-map metric M(r,t), Jacobian factors, and smooth-remainder payloads.\n\nValidation: uv run python scripts/check_docs_freshness.py Entire-Checkpoint: 7ffd824cac89 --- docs/nearfield-template-experiments.md | 104 +++++++++++++++++++++++++ 1 file changed, 104 insertions(+) diff --git a/docs/nearfield-template-experiments.md b/docs/nearfield-template-experiments.md index 070db97..cdfbef7 100644 --- a/docs/nearfield-template-experiments.md +++ b/docs/nearfield-template-experiments.md @@ -84,6 +84,110 @@ The reusable part would be template-space singular model integrals or correction operators. The runtime payload remains map coefficients, Jacobian data, source coefficients, target metadata, and smooth-remainder moment/interpolation data. +## Metric-Field Expansion Tables + +The table-building objective is to separate fixed singular template integrals +from per-chart smooth geometry coefficients. For a self or adjacent chart +interaction, use local source/target coordinates near the singular set and write +`z' = z + delta`. The singular distance model is + +```text +|T(z') - T(z)|^2 = delta^T M(z) delta + higher-order smooth terms, +M(z) = DT(z)^T DT(z). +``` + +Choose a positive-definite reference metric `M0` for one metric bin, chart +family, or local average, and write + +```text +M(z) = M0 + Delta M(z). +``` + +Then the logarithmic singular factor can be expanded as + +```text +log(delta^T M(z) delta) += log(delta^T M0 delta) + + log(1 + (delta^T Delta M(z) delta)/(delta^T M0 delta)). +``` + +If the metric family is binned or normalized so that + +```text +abs((delta^T Delta M(z) delta)/(delta^T M0 delta)) < 1, +``` + +then + +```text +log(delta^T M(z) delta) += log(delta^T M0 delta) + + sum_{k>=1} (-1)^(k+1)/k + ((delta^T Delta M(z) delta)/(delta^T M0 delta))^k. +``` + +Thus the kernel singular part has the schematic expansion + +```text +G_sing(z, z + delta) +~= sum_alpha c_alpha(z) S_alpha(delta; M0), +``` + +where `S_alpha` are fixed singular template functions for the selected reference +metric and `c_alpha(z)` are smooth functions of the entries of `Delta M(z)`. +The expansion can be truncated either by polynomial order in `Delta M`, by +interpolation in metric-entry space, or by a learned/empirical low-rank basis. + +For a bilinear self interaction with template basis functions `psi_i` and +`phi_j`, the singular contribution becomes + +```text +A_ij^sing ~= sum_alpha integral integral + psi_i(z) phi_j(z') c_alpha(z) S_alpha(z'-z; M0) + J(z) J(z') dz' dz. +``` + +Expand the smooth per-chart factor in a template basis `p_beta`: + +```text +c_alpha(z) J(z) J(z') ~= sum_beta q_{alpha beta} p_beta(z, z'). +``` + +Then + +```text +A_ij^sing ~= sum_{alpha,beta} q_{alpha beta} T_{ij alpha beta}, +T_{ij alpha beta} = integral integral + psi_i(z) phi_j(z') p_beta(z, z') S_alpha(z'-z; M0) dz' dz. +``` + +The tensors `T_{ij alpha beta}` are precomputed on fixed template domains. At +runtime, each folded chart only supplies the coefficients `q_{alpha beta}` from +its smooth metric/Jacobian fields and any smooth-remainder representation. + +For the fan map + +```text +T(r,t) = (1-r)V + r C(t), +``` + +the metric entries are explicit: + +```text +T_r = C(t) - V, +T_t = r C'(t), + +M(r,t) = [ |C(t)-V|^2 r (C(t)-V).C'(t) ] + [ r (C(t)-V).C'(t) r^2 |C'(t)|^2 ]. +``` + +For straight, Bezier, or mildly curved polynomial `C(t)`, these entries are +smooth low-parameter functions of `(r,t)`, the seed `V`, and curve coefficients. +This is the mathematical reason precomputed template tables can apply to each +chart piece: the singular table is fixed after choosing `M0`, while observed +folded-decomposition geometry should enter through a small smooth expansion of +`M`, `J`, and the nonsingular remainder. + ## Feasibility Result The answer is qualified yes. Singular or nearly singular folded-piece From 2653068afff0585f099c1d0204322e47e42f40e8 Mon Sep 17 00:00:00 2001 From: xywei Date: Fri, 12 Jun 2026 18:36:41 -0500 Subject: [PATCH 03/15] docs: format near-field derivation math Convert the near-field template derivation from text code fences to display math so GitHub renders the formulas properly. The content is unchanged aside from LaTeX notation cleanup.\n\nValidation: uv run python scripts/check_docs_freshness.py Entire-Checkpoint: c15094b32a03 --- docs/nearfield-template-experiments.md | 178 +++++++++++++------------ 1 file changed, 96 insertions(+), 82 deletions(-) diff --git a/docs/nearfield-template-experiments.md b/docs/nearfield-template-experiments.md index cdfbef7..2fb7a49 100644 --- a/docs/nearfield-template-experiments.md +++ b/docs/nearfield-template-experiments.md @@ -8,17 +8,17 @@ box-code or Volumential work. For one straight-edge folded 2D piece, use the fan chart -```text -T(r, t) = (1-r) V + r C(t), -C(t) = (1-t) P0 + t P1, -0 <= r,t <= 1. -``` +$$ +T(r,t) = (1-r)V + r C(t), +\qquad C(t) = (1-t)P_0 + tP_1, +\qquad 0 \le r,t \le 1. +$$ Its Jacobian is -```text -J(r, t) = r det(C(t)-V, C'(t)) = r det(P0-V, P1-V). -``` +$$ +J(r,t) = r \det(C(t)-V, C'(t)) = r \det(P_0-V, P_1-V). +$$ The factor `r` is the Duffy-style apex degeneracy already used by `cutkit.quadrature.folded2d`. It is not a new singularity for physical @@ -28,32 +28,32 @@ integration; it cancels area at the apex. For the 2D Laplace kernel -```text -G(x, y) = -log(|x-y|)/(2*pi), -``` +$$ +G(x,y) = -\frac{1}{2\pi}\log |x-y|, +$$ the point-target near-field integral is -```text -u(x) = integral_[0,1]^2 G(x, T_e(r,t)) rho(T_e(r,t)) J_e(r,t) dr dt. -``` +$$ +u(x) = \int_{[0,1]^2} G(x,T_e(r,t))\rho(T_e(r,t))J_e(r,t)\,dr\,dt. +$$ This is the cheapest diagnostic path because it isolates source-piece mapping, source density, and target location before introducing target basis functions. For target pieces with a second map `T_f(eta)`, the bilinear template integral is -```text -A_fe = integral_[0,1]^2 integral_[0,1]^2 - psi_f(eta) G(T_f(eta), T_e(xi)) phi_e(xi) - J_f(eta) J_e(xi) dxi deta. -``` +$$ +A_{fe} = \int_{[0,1]^2}\int_{[0,1]^2} +\psi_f(\eta)G(T_f(\eta),T_e(\xi))\phi_e(\xi) +J_f(\eta)J_e(\xi)\,d\xi\,d\eta. +$$ For the self-piece case, this specializes to -```text -I = integral_[0,1]^2 integral_[0,1]^2 - G(T(r,t), T(r',t')) J(r,t) J(r',t') dr dt dr' dt'. -``` +$$ +I = \int_{[0,1]^2}\int_{[0,1]^2} +G(T(r,t),T(r',t'))J(r,t)J(r',t')\,dr\,dt\,dr'\,dt'. +$$ The experiment should classify interactions as: @@ -67,18 +67,22 @@ The experiment should classify interactions as: Near a non-apex diagonal point `z=(r,t)`, write `delta = z' - z`. Then -```text -T(z + delta) - T(z) = DT(z) delta + O(|delta|^2), -|T(z + delta) - T(z)| = sqrt(delta^T M(z) delta) (1 + O(|delta|)), -M(z) = DT(z)^T DT(z). -``` +$$ +T(z+\delta)-T(z) = DT(z)\delta + O(|\delta|^2), +$$ + +$$ +|T(z+\delta)-T(z)| = \sqrt{\delta^T M(z)\delta}\,(1+O(|\delta|)), +\qquad M(z) = DT(z)^TDT(z). +$$ Therefore -```text -G(T(z), T(z+delta)) -= -log(sqrt(delta^T M(z) delta))/(2*pi) + smooth remainder. -``` +$$ +G(T(z),T(z+\delta)) += -\frac{1}{2\pi}\log\sqrt{\delta^T M(z)\delta} ++ \text{smooth remainder}. +$$ The reusable part would be template-space singular model integrals or correction operators. The runtime payload remains map coefficients, Jacobian data, source @@ -91,47 +95,48 @@ from per-chart smooth geometry coefficients. For a self or adjacent chart interaction, use local source/target coordinates near the singular set and write `z' = z + delta`. The singular distance model is -```text -|T(z') - T(z)|^2 = delta^T M(z) delta + higher-order smooth terms, -M(z) = DT(z)^T DT(z). -``` +$$ +|T(z')-T(z)|^2 = \delta^T M(z)\delta + \text{higher-order smooth terms}, +\qquad M(z)=DT(z)^TDT(z). +$$ Choose a positive-definite reference metric `M0` for one metric bin, chart family, or local average, and write -```text -M(z) = M0 + Delta M(z). -``` +$$ +M(z) = M_0 + \Delta M(z). +$$ Then the logarithmic singular factor can be expanded as -```text -log(delta^T M(z) delta) -= log(delta^T M0 delta) - + log(1 + (delta^T Delta M(z) delta)/(delta^T M0 delta)). -``` +$$ +\log(\delta^T M(z)\delta) += \log(\delta^T M_0\delta) ++ \log\left(1+ +\frac{\delta^T\Delta M(z)\delta}{\delta^T M_0\delta}\right). +$$ If the metric family is binned or normalized so that -```text -abs((delta^T Delta M(z) delta)/(delta^T M0 delta)) < 1, -``` +$$ +\left|\frac{\delta^T\Delta M(z)\delta}{\delta^T M_0\delta}\right| < 1, +$$ then -```text -log(delta^T M(z) delta) -= log(delta^T M0 delta) - + sum_{k>=1} (-1)^(k+1)/k - ((delta^T Delta M(z) delta)/(delta^T M0 delta))^k. -``` +$$ +\log(\delta^T M(z)\delta) += \log(\delta^T M_0\delta) ++ \sum_{k\ge 1}\frac{(-1)^{k+1}}{k} +\left(\frac{\delta^T\Delta M(z)\delta}{\delta^T M_0\delta}\right)^k. +$$ Thus the kernel singular part has the schematic expansion -```text -G_sing(z, z + delta) -~= sum_alpha c_alpha(z) S_alpha(delta; M0), -``` +$$ +G_{\mathrm{sing}}(z,z+\delta) +\approx \sum_\alpha c_\alpha(z)S_\alpha(\delta;M_0), +$$ where `S_alpha` are fixed singular template functions for the selected reference metric and `c_alpha(z)` are smooth functions of the entries of `Delta M(z)`. @@ -141,25 +146,29 @@ interpolation in metric-entry space, or by a learned/empirical low-rank basis. For a bilinear self interaction with template basis functions `psi_i` and `phi_j`, the singular contribution becomes -```text -A_ij^sing ~= sum_alpha integral integral - psi_i(z) phi_j(z') c_alpha(z) S_alpha(z'-z; M0) - J(z) J(z') dz' dz. -``` +$$ +A_{ij}^{\mathrm{sing}} \approx \sum_\alpha \int\!\int +\psi_i(z)\phi_j(z')c_\alpha(z)S_\alpha(z'-z;M_0) +J(z)J(z')\,dz'\,dz. +$$ Expand the smooth per-chart factor in a template basis `p_beta`: -```text -c_alpha(z) J(z) J(z') ~= sum_beta q_{alpha beta} p_beta(z, z'). -``` +$$ +c_\alpha(z)J(z)J(z') \approx \sum_\beta q_{\alpha\beta}p_\beta(z,z'). +$$ Then -```text -A_ij^sing ~= sum_{alpha,beta} q_{alpha beta} T_{ij alpha beta}, -T_{ij alpha beta} = integral integral - psi_i(z) phi_j(z') p_beta(z, z') S_alpha(z'-z; M0) dz' dz. -``` +$$ +A_{ij}^{\mathrm{sing}} \approx +\sum_{\alpha,\beta}q_{\alpha\beta}T_{ij\alpha\beta}, +$$ + +$$ +T_{ij\alpha\beta} = \int\!\int +\psi_i(z)\phi_j(z')p_\beta(z,z')S_\alpha(z'-z;M_0)\,dz'\,dz. +$$ The tensors `T_{ij alpha beta}` are precomputed on fixed template domains. At runtime, each folded chart only supplies the coefficients `q_{alpha beta}` from @@ -167,19 +176,24 @@ its smooth metric/Jacobian fields and any smooth-remainder representation. For the fan map -```text -T(r,t) = (1-r)V + r C(t), -``` +$$ +T(r,t) = (1-r)V + rC(t), +$$ the metric entries are explicit: -```text -T_r = C(t) - V, -T_t = r C'(t), +$$ +T_r = C(t)-V, +\qquad T_t = rC'(t), +$$ -M(r,t) = [ |C(t)-V|^2 r (C(t)-V).C'(t) ] - [ r (C(t)-V).C'(t) r^2 |C'(t)|^2 ]. -``` +$$ +M(r,t) = +\begin{bmatrix} +|C(t)-V|^2 & r(C(t)-V)\cdot C'(t) \\ +r(C(t)-V)\cdot C'(t) & r^2|C'(t)|^2 +\end{bmatrix}. +$$ For straight, Bezier, or mildly curved polynomial `C(t)`, these entries are smooth low-parameter functions of `(r,t)`, the seed `V`, and curve coefficients. @@ -247,9 +261,9 @@ uv run python scripts/run_nearfield_template_experiment.py --order 12 The experiment checks the exact scale law for the 2D log kernel. If the fan is scaled by `lambda`, then -```text -I_lambda = lambda^4 (I - log(lambda) A^2/(2*pi)), -``` +$$ +I_\lambda = \lambda^4\left(I - \frac{\log(\lambda)A^2}{2\pi}\right), +$$ where `A` is the signed area of the unscaled fan piece. It also records diagonal remainder samples that shrink as the source/target template coordinates coalesce From 5ffa3027a6eaaf33d5cf92467a7b62814ecbbe97 Mon Sep 17 00:00:00 2001 From: xywei Date: Fri, 12 Jun 2026 18:55:23 -0500 Subject: [PATCH 04/15] fix(evals): handle near-field quadrature edge cases Address Codex review feedback on the near-field template prototype. Reject direct self-interaction quadrature choices with overlapping source/target nodes, report a clear error when point-target quadrature hits a coincident source node, and generalize the log-kernel scale law to density-weighted source and target masses. Also weaken the local metric split documentation so it describes a leading singular model plus bounded/higher-order corrections rather than an unqualified smooth remainder.\n\nValidation: make dev Entire-Checkpoint: 95096c66f419 --- docs/nearfield-template-experiments.md | 24 +++++++--- src/cutkit/evals/__init__.py | 2 + src/cutkit/evals/nearfield_templates.py | 54 ++++++++++++++++++++- tests/evals/test_nearfield_templates.py | 64 +++++++++++++++++++++++++ 4 files changed, 135 insertions(+), 9 deletions(-) diff --git a/docs/nearfield-template-experiments.md b/docs/nearfield-template-experiments.md index 2fb7a49..f1f1a27 100644 --- a/docs/nearfield-template-experiments.md +++ b/docs/nearfield-template-experiments.md @@ -81,12 +81,18 @@ Therefore $$ G(T(z),T(z+\delta)) = -\frac{1}{2\pi}\log\sqrt{\delta^T M(z)\delta} -+ \text{smooth remainder}. ++ \text{lower-order correction}. $$ The reusable part would be template-space singular model integrals or correction operators. The runtime payload remains map coefficients, Jacobian data, source -coefficients, target metadata, and smooth-remainder moment/interpolation data. +coefficients, target metadata, and correction moment/interpolation data. For the +straight fan map, `T(r,t)` is bilinear in template variables, so subtracting only +the local metric model leaves an `O(|delta|)` correction whose first derivative +at the diagonal can depend on approach direction. A smooth-remainder expansion +should therefore either include higher-order distance terms in the singular model +or treat this first prototype as a leading singular split plus a bounded local +correction. ## Metric-Field Expansion Tables @@ -96,7 +102,7 @@ interaction, use local source/target coordinates near the singular set and write `z' = z + delta`. The singular distance model is $$ -|T(z')-T(z)|^2 = \delta^T M(z)\delta + \text{higher-order smooth terms}, +|T(z')-T(z)|^2 = \delta^T M(z)\delta + \text{higher-order chart terms}, \qquad M(z)=DT(z)^TDT(z). $$ @@ -142,6 +148,9 @@ where `S_alpha` are fixed singular template functions for the selected reference metric and `c_alpha(z)` are smooth functions of the entries of `Delta M(z)`. The expansion can be truncated either by polynomial order in `Delta M`, by interpolation in metric-entry space, or by a learned/empirical low-rank basis. +For full smooth-remainder tables near a chart diagonal, the singular model should +also include enough higher-order chart-distance terms to remove direction- +dependent local corrections. For a bilinear self interaction with template basis functions `psi_i` and `phi_j`, the singular contribution becomes @@ -172,7 +181,8 @@ $$ The tensors `T_{ij alpha beta}` are precomputed on fixed template domains. At runtime, each folded chart only supplies the coefficients `q_{alpha beta}` from -its smooth metric/Jacobian fields and any smooth-remainder representation. +its smooth metric/Jacobian fields and any higher-order correction or remainder +representation. For the fan map @@ -200,7 +210,7 @@ smooth low-parameter functions of `(r,t)`, the seed `V`, and curve coefficients. This is the mathematical reason precomputed template tables can apply to each chart piece: the singular table is fixed after choosing `M0`, while observed folded-decomposition geometry should enter through a small smooth expansion of -`M`, `J`, and the nonsingular remainder. +`M`, `J`, and the higher-order correction data. ## Feasibility Result @@ -215,7 +225,7 @@ because the folded maps have smooth, low-parameter structure away from collapsed seed faces, the metric field `M(z)` should vary smoothly over the template and across the geometry families produced by folded decomposition. A functional expansion in the metric data, for example moments, interpolation coefficients, -or a low-rank basis in `M` and the smooth remainder, may cover enough practical +or a low-rank basis in `M` and the correction terms, may cover enough practical cut-piece cases to make reusable near-field corrections worthwhile. The open numerical question is therefore whether the smooth dependence on `M`, @@ -229,7 +239,7 @@ For each geometry and interaction case, record: - reference value; - ordinary pulled-forward tensor-product quadrature error; - template singular-correction error; -- smooth-remainder approximation error; +- higher-order correction/remainder approximation error; - dependence on quadrature order; - dependence on seed location and curve coefficients; - size, rank, and smoothness of geometry-parameterized correction data; diff --git a/src/cutkit/evals/__init__.py b/src/cutkit/evals/__init__.py index 245c6d2..a1f1557 100644 --- a/src/cutkit/evals/__init__.py +++ b/src/cutkit/evals/__init__.py @@ -85,6 +85,7 @@ point_target_laplace_potential, run_nearfield_template_experiment, self_interaction_laplace, + template_density_mass, ) __all__ = [ @@ -160,6 +161,7 @@ "point_target_laplace_potential", "run_nearfield_template_experiment", "self_interaction_laplace", + "template_density_mass", "FormDslParityBenchmark", "FormDslParityRow", "FormDslMultipatchStressBenchmark", diff --git a/src/cutkit/evals/nearfield_templates.py b/src/cutkit/evals/nearfield_templates.py index d28e5d7..b4ec56a 100644 --- a/src/cutkit/evals/nearfield_templates.py +++ b/src/cutkit/evals/nearfield_templates.py @@ -10,6 +10,7 @@ Point2D = tuple[float, float] TemplateDensity2D = Callable[[float, float], float] +_COINCIDENT_TOL = 1.0e-14 def unit_template_density(_r: float, _t: float) -> float: @@ -34,6 +35,19 @@ def _scale_point(point: Point2D, scale: float) -> Point2D: return (scale * point[0], scale * point[1]) +def _rules_share_nodes( + lhs: tuple[float, ...], + rhs: tuple[float, ...], + *, + atol: float = _COINCIDENT_TOL, +) -> bool: + for left in lhs: + for right in rhs: + if abs(left - right) <= atol: + return True + return False + + @dataclass(frozen=True) class FanTemplateMap2D: """Straight-edge folded fan chart ``T(r, t) = (1-r)V + r C(t)``.""" @@ -183,6 +197,11 @@ def self_interaction_laplace( target_nodes, target_weights = gauss_legendre_01(order) source_nodes, source_weights = gauss_legendre_01(source_order) + if _rules_share_nodes(target_nodes, source_nodes): + raise ValueError( + "source_order quadrature nodes must not overlap target order nodes " + "for direct self-interaction reference quadrature" + ) total = 0.0 for target_r, target_wr in zip(target_nodes, target_weights, strict=True): @@ -213,6 +232,28 @@ def self_interaction_laplace( return total +def template_density_mass( + fan: FanTemplateMap2D, + *, + order: int, + density: TemplateDensity2D | None = None, +) -> float: + """Return ``integral density(r,t) * J(r,t) dr dt`` on one fan chart.""" + + if order < 1: + raise ValueError("order must be positive") + if density is None: + density = unit_template_density + + nodes, weights = gauss_legendre_01(order) + total = 0.0 + for r, wr in zip(nodes, weights, strict=True): + jacobian = fan.signed_jacobian(r) + for t, wt in zip(nodes, weights, strict=True): + total += density(r, t) * jacobian * wr * wt + return total + + def point_target_laplace_potential( fan: FanTemplateMap2D, target: Point2D, @@ -233,6 +274,11 @@ def point_target_laplace_potential( jacobian = fan.signed_jacobian(r) for t, wt in zip(nodes, weights, strict=True): source = fan.point(r, t) + if hypot(target[0] - source[0], target[1] - source[1]) <= _COINCIDENT_TOL: + raise ValueError( + "point target coincides with a source quadrature node; " + "use a singular point-target reference rule" + ) total += ( laplace_log_kernel(target, source) * source_density(r, t) @@ -245,15 +291,19 @@ def point_target_laplace_potential( def expected_scaled_laplace_self_interaction( base_interaction: float, - signed_area: float, + source_mass: float, scale_factor: float, + *, + target_mass: float | None = None, ) -> float: """Return the exact 2D log-kernel scale law for a scaled fan piece.""" if scale_factor <= 0.0: raise ValueError("scale_factor must be positive") + if target_mass is None: + target_mass = source_mass return scale_factor**4 * ( - base_interaction - log(scale_factor) * signed_area * signed_area / (2.0 * pi) + base_interaction - log(scale_factor) * target_mass * source_mass / (2.0 * pi) ) diff --git a/tests/evals/test_nearfield_templates.py b/tests/evals/test_nearfield_templates.py index 23cb351..126a9e2 100644 --- a/tests/evals/test_nearfield_templates.py +++ b/tests/evals/test_nearfield_templates.py @@ -9,7 +9,9 @@ point_target_laplace_potential, run_nearfield_template_experiment, self_interaction_laplace, + template_density_mass, ) +from cutkit.quadrature import gauss_legendre_01 def test_fan_map_area_and_metric_are_consistent() -> None: @@ -48,6 +50,55 @@ def test_self_interaction_obeys_log_kernel_scale_law() -> None: assert scaled == pytest.approx(expected, abs=1.0e-13) +def test_self_interaction_rejects_overlapping_template_nodes() -> None: + fan = FanTemplateMap2D( + vertex=(0.0, 0.0), + edge_start=(1.0, 0.0), + edge_end=(0.35, 0.9), + ) + + with pytest.raises(ValueError, match="must not overlap"): + self_interaction_laplace(fan, order=3, source_order=5) + + +def test_scaled_self_interaction_uses_density_weighted_masses() -> None: + fan = FanTemplateMap2D( + vertex=(0.0, 0.0), + edge_start=(1.0, 0.0), + edge_end=(0.35, 0.9), + ) + + def source_density(r: float, t: float) -> float: + return 1.0 + r + + def target_density(r: float, t: float) -> float: + return 1.0 - 0.25 * t + + scale_factor = 1.4 + value = self_interaction_laplace( + fan, + order=4, + source_order=5, + source_density=source_density, + target_density=target_density, + ) + scaled = self_interaction_laplace( + fan.scaled(scale_factor), + order=4, + source_order=5, + source_density=source_density, + target_density=target_density, + ) + expected = expected_scaled_laplace_self_interaction( + value, + template_density_mass(fan, order=4, density=source_density), + scale_factor, + target_mass=template_density_mass(fan, order=4, density=target_density), + ) + + assert scaled == pytest.approx(expected, abs=1.0e-13) + + def test_metric_singular_remainder_shrinks_toward_diagonal() -> None: fan = FanTemplateMap2D( vertex=(0.0, 0.0), @@ -94,6 +145,19 @@ def density(r: float, t: float) -> float: assert abs(fine - reference) < abs(coarse - reference) +def test_point_target_path_rejects_coincident_source_node() -> None: + fan = FanTemplateMap2D( + vertex=(0.0, 0.0), + edge_start=(1.0, 0.0), + edge_end=(0.35, 0.9), + ) + nodes, _weights = gauss_legendre_01(3) + target = fan.point(nodes[1], nodes[1]) + + with pytest.raises(ValueError, match="coincides with a source quadrature node"): + point_target_laplace_potential(fan, target, order=3) + + def test_baseline_experiment_records_feasibility_checks() -> None: result = run_nearfield_template_experiment(order=5, scale_factor=1.5) From 38560bbfc3d14dc707119593268798cfbe1baf21 Mon Sep 17 00:00:00 2001 From: xywei Date: Fri, 12 Jun 2026 18:58:32 -0500 Subject: [PATCH 05/15] docs: make Bezier fan chart the primary derivation Update the near-field template derivation to start from a CAD-realistic Bezier trim fan chart instead of a straight-edge specialization. Straight fan pieces are now described as smoke-test fixtures, while the metric-field expansion is framed in terms of Bezier control points and seed locations.\n\nValidation: uv run python scripts/check_docs_freshness.py Entire-Checkpoint: c9e0b70886a7 --- docs/nearfield-template-experiments.md | 69 +++++++++++++++----------- 1 file changed, 41 insertions(+), 28 deletions(-) diff --git a/docs/nearfield-template-experiments.md b/docs/nearfield-template-experiments.md index f1f1a27..a182b3f 100644 --- a/docs/nearfield-template-experiments.md +++ b/docs/nearfield-template-experiments.md @@ -4,25 +4,36 @@ This note records the first mathematical check for issue 61: whether folded decomposition charts can support a reusable near-field template story for later box-code or Volumential work. -## Fan Map +## Bezier Fan Chart -For one straight-edge folded 2D piece, use the fan chart +For the mathematical model, start with the chart type expected from CAD-backed +folded decomposition: a fan from a seed point `V` to a smooth Bezier trim curve +`C(t)`. A degree-`p` Bezier edge is $$ -T(r,t) = (1-r)V + r C(t), -\qquad C(t) = (1-t)P_0 + tP_1, -\qquad 0 \le r,t \le 1. +C(t) = \sum_{a=0}^p B_a^p(t)P_a, +\qquad B_a^p(t)=\binom{p}{a}(1-t)^{p-a}t^a, +\qquad 0\le t\le 1. +$$ + +The folded chart is + +$$ +T(r,t) = (1-r)V + rC(t), +\qquad 0\le r,t\le 1. $$ Its Jacobian is $$ -J(r,t) = r \det(C(t)-V, C'(t)) = r \det(P_0-V, P_1-V). +J(r,t) = r\det(C(t)-V,C'(t)). $$ The factor `r` is the Duffy-style apex degeneracy already used by `cutkit.quadrature.folded2d`. It is not a new singularity for physical -integration; it cancels area at the apex. +integration; it cancels area at the apex. The straight-edge case is only the +degree-1 specialization and should be treated as a smoke-test fixture, not as +the primary derivation. ## Mapped Kernel @@ -86,13 +97,13 @@ $$ The reusable part would be template-space singular model integrals or correction operators. The runtime payload remains map coefficients, Jacobian data, source -coefficients, target metadata, and correction moment/interpolation data. For the -straight fan map, `T(r,t)` is bilinear in template variables, so subtracting only -the local metric model leaves an `O(|delta|)` correction whose first derivative -at the diagonal can depend on approach direction. A smooth-remainder expansion -should therefore either include higher-order distance terms in the singular model -or treat this first prototype as a leading singular split plus a bounded local -correction. +coefficients, target metadata, and correction moment/interpolation data. For a +Bezier fan map, `T(r,t)` is polynomial in `(r,t)` and contains mixed higher-order +terms from `rC(t)`. Subtracting only the local metric model leaves an +`O(|delta|)` correction whose first derivative at the diagonal can depend on +approach direction. A smooth-remainder expansion should therefore either include +higher-order distance terms in the singular model or treat this first prototype +as a leading singular split plus a bounded local correction. ## Metric-Field Expansion Tables @@ -184,7 +195,7 @@ runtime, each folded chart only supplies the coefficients `q_{alpha beta}` from its smooth metric/Jacobian fields and any higher-order correction or remainder representation. -For the fan map +For the Bezier fan map $$ T(r,t) = (1-r)V + rC(t), @@ -205,12 +216,12 @@ r(C(t)-V)\cdot C'(t) & r^2|C'(t)|^2 \end{bmatrix}. $$ -For straight, Bezier, or mildly curved polynomial `C(t)`, these entries are -smooth low-parameter functions of `(r,t)`, the seed `V`, and curve coefficients. -This is the mathematical reason precomputed template tables can apply to each -chart piece: the singular table is fixed after choosing `M0`, while observed -folded-decomposition geometry should enter through a small smooth expansion of -`M`, `J`, and the higher-order correction data. +For Bezier or piecewise-Bezier CAD trims, these entries are smooth +low-parameter functions of `(r,t)`, the seed `V`, and the Bezier control points +`P_a`. This is the mathematical reason precomputed template tables can apply to +each chart piece: the singular table is fixed after choosing `M0`, while +observed folded-decomposition geometry should enter through a small smooth +expansion of `M`, `J`, and the higher-order correction data. ## Feasibility Result @@ -246,16 +257,18 @@ For each geometry and interaction case, record: - how many metric-field expansion modes are needed for the observed folded decomposition cases. -The first geometry family is CAD-independent: straight fan pieces, then -quadratic Bezier arcs and mildly curved cubics, with seed locations chosen to -produce positive, negative, and folded orientations. OpenCascade should only be -used later as a stress-test backend. +The first geometry family should be CAD-independent but CAD-realistic: quadratic +and cubic Bezier fan charts with seed locations chosen to produce positive, +negative, and folded orientations. Straight fan pieces are useful only as +low-level smoke tests. OpenCascade should be used later as a source of +stress-test Bezier/NURBS-derived fixtures, not as a dependency of the numerical +question. ## Prototype `cutkit.evals.nearfield_templates` implements: -- `FanTemplateMap2D` for the analytic fan map. +- `FanTemplateMap2D` for the current straight-edge smoke-test fan map. - `point_target_laplace_potential(...)` for point-target source integrals with polynomial template densities. - `self_interaction_laplace(...)` for a direct high-order mapped self integral. @@ -268,8 +281,8 @@ The script wrapper is: uv run python scripts/run_nearfield_template_experiment.py --order 12 ``` -The experiment checks the exact scale law for the 2D log kernel. If the fan is -scaled by `lambda`, then +The current smoke-test experiment checks the exact scale law for the 2D log +kernel. If the fan is scaled by `lambda`, then $$ I_\lambda = \lambda^4\left(I - \frac{\log(\lambda)A^2}{2\pi}\right), From db19ee4e0dfe575ed2aa22c00170fcae7c90e671 Mon Sep 17 00:00:00 2001 From: xywei Date: Fri, 12 Jun 2026 19:25:40 -0500 Subject: [PATCH 06/15] fix(evals): make near-field prototype point-target only Remove the folded-target and bilinear self-interaction framing from the near-field template experiment. The docs now describe source folded pieces evaluated at physical target nodes in containing and neighboring boxes, and the experiment API now reports point-target scale-law quantities only.\n\nValidation: make dev; uv run --extra dev python -m pytest tests/evals/test_nearfield_templates.py tests/test_script_wrappers.py; uv run python scripts/check_docs_freshness.py Entire-Checkpoint: 6cb8d3064a3b --- ...2-boxcode-nearfield-template-experiment.md | 12 +- docs/nearfield-template-experiments.md | 140 ++++++++++-------- scripts/run_nearfield_template_experiment.py | 11 +- src/cutkit/evals/__init__.py | 6 +- src/cutkit/evals/nearfield_templates.py | 132 ++++------------- tests/evals/test_nearfield_templates.py | 52 +++---- 6 files changed, 143 insertions(+), 210 deletions(-) diff --git a/docs/exec-plans/completed/20260612-boxcode-nearfield-template-experiment.md b/docs/exec-plans/completed/20260612-boxcode-nearfield-template-experiment.md index 5c85bee..4b07f66 100644 --- a/docs/exec-plans/completed/20260612-boxcode-nearfield-template-experiment.md +++ b/docs/exec-plans/completed/20260612-boxcode-nearfield-template-experiment.md @@ -8,7 +8,7 @@ template corrections for box-code and Volumential-style workflows. ## Scope - Derive the mapped-kernel form for a 2D fan folded piece. -- Include both point-target and target-piece bilinear near-field forms. +- Focus the near-field form on point targets in physical space. - Identify singular factors that can be separated from smooth geometry payloads. - Frame the reuse hypothesis as a functional expansion in smooth metric fields `M(z)`, not just as raw runtime geometry data. @@ -24,9 +24,9 @@ template corrections for box-code and Volumential-style workflows. ## Acceptance Criteria - The derivation states the mapped self-interaction form and its split. -- The derivation records the point-target form, bilinear form, interaction - taxonomy, and measurement criteria from issue 61. -- The prototype records a self or adjacent template experiment with a direct +- The derivation records the point-target form, target/source geometry taxonomy, + and measurement criteria from issue 61. +- The prototype records a point-target template experiment with a direct high-order reference value. - Tests check the experiment against an analytic or independently computed invariant. @@ -53,8 +53,8 @@ template corrections for box-code and Volumential-style workflows. - The practical reuse hypothesis is that a functional expansion in `M(z)` and smooth-remainder data covers the folded-decomposition cases with few modes. - The first prototype covers analytic straight fan maps, point-target - integrals, bilinear self interactions, a log-kernel scale-law check, and a - diagonal metric-remainder sample. + integrals, a log-kernel scale-law check, and a diagonal metric-remainder + sample. ## Validation diff --git a/docs/nearfield-template-experiments.md b/docs/nearfield-template-experiments.md index a182b3f..5ef90e3 100644 --- a/docs/nearfield-template-experiments.md +++ b/docs/nearfield-template-experiments.md @@ -49,72 +49,87 @@ $$ u(x) = \int_{[0,1]^2} G(x,T_e(r,t))\rho(T_e(r,t))J_e(r,t)\,dr\,dt. $$ -This is the cheapest diagnostic path because it isolates source-piece mapping, -source density, and target location before introducing target basis functions. -For target pieces with a second map `T_f(eta)`, the bilinear template integral is +This is the primary target model for the box-code experiment. The source is a +folded piece, but targets are physical-space discretization nodes in the source +box and neighboring near-field boxes. There is no target folded piece in the +first near-field experiment. + +For reusable tables, express the target relative to the same physical box or +chart payload. If `x0` is a chart/box reference point and `H` is a box-scale map, +write $$ -A_{fe} = \int_{[0,1]^2}\int_{[0,1]^2} -\psi_f(\eta)G(T_f(\eta),T_e(\xi))\phi_e(\xi) -J_f(\eta)J_e(\xi)\,d\xi\,d\eta. +x = x_0 + H\tau, $$ -For the self-piece case, this specializes to +where `tau` is the template-space location of a target node in the containing or +neighboring box. The source-to-point integral is then a fixed-domain integral in +the source chart coordinates and the target offset parameter `tau`: $$ -I = \int_{[0,1]^2}\int_{[0,1]^2} -G(T(r,t),T(r',t'))J(r,t)J(r',t')\,dr\,dt\,dr'\,dt'. +u(\tau) = \int_{[0,1]^2}G(x_0+H\tau,T_e(\xi))\rho(T_e(\xi))J_e(\xi)\,d\xi. $$ -The experiment should classify interactions as: +The experiment should classify target/source geometry as: -- self piece, `e = f`; -- edge-adjacent pieces sharing a boundary segment; -- vertex or seed-adjacent pieces; -- near but disjoint pieces; -- well-separated pieces, where ordinary tensor-product quadrature should already work. +- target node on or very near the source folded piece; +- target node in the source piece's containing box; +- target node in an edge-neighboring box; +- target node in a vertex-neighboring box; +- well-separated target node, where ordinary tensor-product quadrature should already work. -## Singular Split +## Point-Target Singular Split -Near a non-apex diagonal point `z=(r,t)`, write `delta = z' - z`. Then +For a target point near the source chart, let `z_x` be a closest or projected +template coordinate with `T(z_x)` near `x`, and write `xi = z_x + delta`. Then $$ -T(z+\delta)-T(z) = DT(z)\delta + O(|\delta|^2), +T(z_x+\delta)-T(z_x) = DT(z_x)\delta + O(|\delta|^2), $$ $$ -|T(z+\delta)-T(z)| = \sqrt{\delta^T M(z)\delta}\,(1+O(|\delta|)), -\qquad M(z) = DT(z)^TDT(z). +|T(z_x+\delta)-T(z_x)| = +\sqrt{\delta^T M(z_x)\delta}\,(1+O(|\delta|)), +\qquad M(z_x) = DT(z_x)^TDT(z_x). $$ -Therefore +For an on-surface or asymptotically close target, this gives the leading model $$ -G(T(z),T(z+\delta)) -= -\frac{1}{2\pi}\log\sqrt{\delta^T M(z)\delta} +G(x,T(z_x+\delta)) += -\frac{1}{2\pi}\log\sqrt{\delta^T M(z_x)\delta} + \text{lower-order correction}. $$ +For an off-surface target, include the normal/offset residual `d = x-T(z_x)`: + +$$ +|x-T(z_x+\delta)|^2 += |d-DT(z_x)\delta|^2 + \text{higher-order chart terms}. +$$ + The reusable part would be template-space singular model integrals or correction operators. The runtime payload remains map coefficients, Jacobian data, source -coefficients, target metadata, and correction moment/interpolation data. For a -Bezier fan map, `T(r,t)` is polynomial in `(r,t)` and contains mixed higher-order -terms from `rC(t)`. Subtracting only the local metric model leaves an -`O(|delta|)` correction whose first derivative at the diagonal can depend on -approach direction. A smooth-remainder expansion should therefore either include -higher-order distance terms in the singular model or treat this first prototype -as a leading singular split plus a bounded local correction. +coefficients, target-node offsets, and correction moment/interpolation data. For +a Bezier fan map, `T(r,t)` is polynomial in `(r,t)` and contains mixed +higher-order terms from `rC(t)`. Subtracting only the local metric model leaves +an `O(|delta|)` correction whose first derivative at the singular point can +depend on approach direction. A smooth-remainder expansion should therefore +either include higher-order distance terms in the singular model or treat this +first prototype as a leading singular split plus a bounded local correction. ## Metric-Field Expansion Tables The table-building objective is to separate fixed singular template integrals -from per-chart smooth geometry coefficients. For a self or adjacent chart -interaction, use local source/target coordinates near the singular set and write -`z' = z + delta`. The singular distance model is +from per-chart smooth geometry and per-target offset coefficients. For a target +node represented by offset coordinates `tau`, use local source coordinates near +the closest chart point and write `xi = z_x + delta`. The on-surface singular +distance model is $$ -|T(z')-T(z)|^2 = \delta^T M(z)\delta + \text{higher-order chart terms}, -\qquad M(z)=DT(z)^TDT(z). +|T(z_x+\delta)-T(z_x)|^2 += \delta^T M(z_x)\delta + \text{higher-order chart terms}, +\qquad M(z_x)=DT(z_x)^TDT(z_x). $$ Choose a positive-definite reference metric `M0` for one metric bin, chart @@ -151,49 +166,50 @@ $$ Thus the kernel singular part has the schematic expansion $$ -G_{\mathrm{sing}}(z,z+\delta) -\approx \sum_\alpha c_\alpha(z)S_\alpha(\delta;M_0), +G_{\mathrm{sing}}(x,T(z_x+\delta)) +\approx \sum_\alpha c_\alpha(z_x,\tau)S_\alpha(\delta;M_0), $$ where `S_alpha` are fixed singular template functions for the selected reference -metric and `c_alpha(z)` are smooth functions of the entries of `Delta M(z)`. -The expansion can be truncated either by polynomial order in `Delta M`, by -interpolation in metric-entry space, or by a learned/empirical low-rank basis. +metric and `c_alpha(z_x,tau)` are smooth functions of the entries of +`Delta M(z_x)` and the target offset data. The expansion can be truncated either +by polynomial order in `Delta M`, by interpolation in metric/offset space, or by +a learned/empirical low-rank basis. For full smooth-remainder tables near a chart diagonal, the singular model should also include enough higher-order chart-distance terms to remove direction- dependent local corrections. -For a bilinear self interaction with template basis functions `psi_i` and -`phi_j`, the singular contribution becomes +For one target node and one source density expansion mode `rho_j`, the singular +contribution becomes $$ -A_{ij}^{\mathrm{sing}} \approx \sum_\alpha \int\!\int -\psi_i(z)\phi_j(z')c_\alpha(z)S_\alpha(z'-z;M_0) -J(z)J(z')\,dz'\,dz. +u_j^{\mathrm{sing}}(\tau) \approx \sum_\alpha \int_{[0,1]^2} +\rho_j(\xi)c_\alpha(z_x,\tau)S_\alpha(\xi-z_x;M_0)J(\xi)\,d\xi. $$ -Expand the smooth per-chart factor in a template basis `p_beta`: +Expand the smooth per-chart/per-target factor in a template basis `p_beta`: $$ -c_\alpha(z)J(z)J(z') \approx \sum_\beta q_{\alpha\beta}p_\beta(z,z'). +c_\alpha(z_x,\tau)J(\xi) \approx \sum_\beta q_{\alpha\beta}(\tau)p_\beta(\xi). $$ Then $$ -A_{ij}^{\mathrm{sing}} \approx -\sum_{\alpha,\beta}q_{\alpha\beta}T_{ij\alpha\beta}, +u_j^{\mathrm{sing}}(\tau) \approx +\sum_{\alpha,\beta}q_{\alpha\beta}(\tau)T_{j\alpha\beta}, $$ $$ -T_{ij\alpha\beta} = \int\!\int -\psi_i(z)\phi_j(z')p_\beta(z,z')S_\alpha(z'-z;M_0)\,dz'\,dz. +T_{j\alpha\beta} = \int_{[0,1]^2} +\rho_j(\xi)p_\beta(\xi)S_\alpha(\xi-z_x;M_0)\,d\xi. $$ -The tensors `T_{ij alpha beta}` are precomputed on fixed template domains. At -runtime, each folded chart only supplies the coefficients `q_{alpha beta}` from -its smooth metric/Jacobian fields and any higher-order correction or remainder -representation. +The tensors `T_{j alpha beta}` are precomputed on fixed source template domains, +or tabulated over a small target-offset grid. At runtime, each folded chart and +target node only supply the coefficients `q_{alpha beta}(tau)` from smooth +metric/Jacobian fields, target offset data, and any higher-order correction or +remainder representation. For the Bezier fan map @@ -271,7 +287,6 @@ question. - `FanTemplateMap2D` for the current straight-edge smoke-test fan map. - `point_target_laplace_potential(...)` for point-target source integrals with polynomial template densities. -- `self_interaction_laplace(...)` for a direct high-order mapped self integral. - `diagonal_remainder_sample(...)` for the metric singular split. - `run_nearfield_template_experiment(...)` for the baseline feasibility check. @@ -282,15 +297,18 @@ uv run python scripts/run_nearfield_template_experiment.py --order 12 ``` The current smoke-test experiment checks the exact scale law for the 2D log -kernel. If the fan is scaled by `lambda`, then +kernel. If both the source fan and physical target point are scaled by `lambda`, +then $$ -I_\lambda = \lambda^4\left(I - \frac{\log(\lambda)A^2}{2\pi}\right), +u_\lambda(\lambda x) += \lambda^2\left(u(x) - \frac{\log(\lambda)m_\rho}{2\pi}\right), $$ -where `A` is the signed area of the unscaled fan piece. It also records diagonal -remainder samples that shrink as the source/target template coordinates coalesce -and a point-target low-order-vs-reference error. +where `m_rho` is the density-weighted signed source mass on the unscaled fan +piece. It also records diagonal remainder samples that shrink as the source and +target template coordinates coalesce and a point-target low-order-vs-reference +error. ## Next Decision diff --git a/scripts/run_nearfield_template_experiment.py b/scripts/run_nearfield_template_experiment.py index 1e4442f..4e5991f 100755 --- a/scripts/run_nearfield_template_experiment.py +++ b/scripts/run_nearfield_template_experiment.py @@ -35,17 +35,20 @@ def main() -> int: print("field | value") print("--- | ---") print(f"order | {result.order}") - print(f"source_order | {result.source_order}") print(f"near_point_target | {result.near_point_target}") print(f"point_target_reference | {result.point_target_reference:.16e}") print(f"point_target_low_order | {result.point_target_low_order:.16e}") print(f"point_target_abs_error | {result.point_target_abs_error:.16e}") print(f"signed_area | {result.signed_area:.16e}") - print(f"self_interaction | {result.self_interaction:.16e}") + print(f"density_mass | {result.density_mass:.16e}") print(f"scale_factor | {result.scale_factor:.16e}") - print(f"scaled_self_interaction | {result.scaled_self_interaction:.16e}") + print(f"scaled_near_point_target | {result.scaled_near_point_target}") print( - f"expected_scaled_self_interaction | {result.expected_scaled_self_interaction:.16e}" + f"scaled_point_target_potential | {result.scaled_point_target_potential:.16e}" + ) + print( + "expected_scaled_point_target_potential | " + f"{result.expected_scaled_point_target_potential:.16e}" ) print(f"scaled_abs_error | {result.scaled_abs_error:.16e}") print("") diff --git a/src/cutkit/evals/__init__.py b/src/cutkit/evals/__init__.py index a1f1557..1b9956b 100644 --- a/src/cutkit/evals/__init__.py +++ b/src/cutkit/evals/__init__.py @@ -79,12 +79,11 @@ FanTemplateMap2D, NearfieldTemplateExperiment, diagonal_remainder_sample, - expected_scaled_laplace_self_interaction, + expected_scaled_laplace_point_potential, laplace_log_kernel, metric_model_distance, point_target_laplace_potential, run_nearfield_template_experiment, - self_interaction_laplace, template_density_mass, ) @@ -155,12 +154,11 @@ "FanTemplateMap2D", "NearfieldTemplateExperiment", "diagonal_remainder_sample", - "expected_scaled_laplace_self_interaction", + "expected_scaled_laplace_point_potential", "laplace_log_kernel", "metric_model_distance", "point_target_laplace_potential", "run_nearfield_template_experiment", - "self_interaction_laplace", "template_density_mass", "FormDslParityBenchmark", "FormDslParityRow", diff --git a/src/cutkit/evals/nearfield_templates.py b/src/cutkit/evals/nearfield_templates.py index b4ec56a..0fdd49f 100644 --- a/src/cutkit/evals/nearfield_templates.py +++ b/src/cutkit/evals/nearfield_templates.py @@ -35,19 +35,6 @@ def _scale_point(point: Point2D, scale: float) -> Point2D: return (scale * point[0], scale * point[1]) -def _rules_share_nodes( - lhs: tuple[float, ...], - rhs: tuple[float, ...], - *, - atol: float = _COINCIDENT_TOL, -) -> bool: - for left in lhs: - for right in rhs: - if abs(left - right) <= atol: - return True - return False - - @dataclass(frozen=True) class FanTemplateMap2D: """Straight-edge folded fan chart ``T(r, t) = (1-r)V + r C(t)``.""" @@ -126,21 +113,21 @@ class DiagonalRemainderSample: @dataclass(frozen=True) class NearfieldTemplateExperiment: - """Summary of one fan-chart self-interaction experiment.""" + """Summary of one fan-chart point-target experiment.""" fan: FanTemplateMap2D order: int - source_order: int near_point_target: Point2D point_target_reference: float point_target_low_order: float point_target_abs_error: float - self_interaction: float + scaled_near_point_target: Point2D + scaled_point_target_potential: float + expected_scaled_point_target_potential: float + scaled_abs_error: float signed_area: float + density_mass: float scale_factor: float - scaled_self_interaction: float - expected_scaled_self_interaction: float - scaled_abs_error: float diagonal_remainders: tuple[DiagonalRemainderSample, ...] @@ -174,64 +161,6 @@ def metric_model_distance( return sqrt(distance_squared) -def self_interaction_laplace( - fan: FanTemplateMap2D, - *, - order: int, - source_order: int | None = None, - source_density: TemplateDensity2D | None = None, - target_density: TemplateDensity2D | None = None, -) -> float: - """Approximate the mapped self interaction on one fan chart.""" - - if order < 1: - raise ValueError("order must be positive") - if source_order is None: - source_order = order + 1 - if source_order < 1: - raise ValueError("source_order must be positive") - if source_density is None: - source_density = unit_template_density - if target_density is None: - target_density = unit_template_density - - target_nodes, target_weights = gauss_legendre_01(order) - source_nodes, source_weights = gauss_legendre_01(source_order) - if _rules_share_nodes(target_nodes, source_nodes): - raise ValueError( - "source_order quadrature nodes must not overlap target order nodes " - "for direct self-interaction reference quadrature" - ) - total = 0.0 - - for target_r, target_wr in zip(target_nodes, target_weights, strict=True): - target_j = fan.signed_jacobian(target_r) - for target_t, target_wt in zip(target_nodes, target_weights, strict=True): - target = fan.point(target_r, target_t) - target_weight = ( - target_density(target_r, target_t) * target_j * target_wr * target_wt - ) - for source_r, source_wr in zip(source_nodes, source_weights, strict=True): - source_j = fan.signed_jacobian(source_r) - for source_t, source_wt in zip( - source_nodes, source_weights, strict=True - ): - source = fan.point(source_r, source_t) - source_weight = ( - source_density(source_r, source_t) - * source_j - * source_wr - * source_wt - ) - total += ( - laplace_log_kernel(target, source) - * target_weight - * source_weight - ) - - return total - - def template_density_mass( fan: FanTemplateMap2D, *, @@ -289,21 +218,17 @@ def point_target_laplace_potential( return total -def expected_scaled_laplace_self_interaction( - base_interaction: float, +def expected_scaled_laplace_point_potential( + base_potential: float, source_mass: float, scale_factor: float, - *, - target_mass: float | None = None, ) -> float: - """Return the exact 2D log-kernel scale law for a scaled fan piece.""" + """Return the exact 2D log-kernel scale law for a scaled source and target.""" if scale_factor <= 0.0: raise ValueError("scale_factor must be positive") - if target_mass is None: - target_mass = source_mass - return scale_factor**4 * ( - base_interaction - log(scale_factor) * target_mass * source_mass / (2.0 * pi) + return scale_factor**2 * ( + base_potential - log(scale_factor) * source_mass / (2.0 * pi) ) @@ -355,12 +280,6 @@ def run_nearfield_template_experiment( edge_start=(1.0, 0.0), edge_end=(0.35, 0.9), ) - source_order = order + 1 - self_value = self_interaction_laplace( - fan, - order=order, - source_order=source_order, - ) near_point_target = (0.47, 0.42) def point_density(r: float, t: float) -> float: @@ -378,14 +297,21 @@ def point_density(r: float, t: float) -> float: order=max(order // 2, 2), source_density=point_density, ) - scaled_value = self_interaction_laplace( + scaled_target = _scale_point(near_point_target, scale_factor) + scaled_value = point_target_laplace_potential( fan.scaled(scale_factor), - order=order, - source_order=source_order, + scaled_target, + order=max(order + 8, 16), + source_density=point_density, + ) + density_mass = template_density_mass( + fan, + order=max(order + 8, 16), + density=point_density, ) - expected_scaled = expected_scaled_laplace_self_interaction( - self_value, - fan.signed_area, + expected_scaled = expected_scaled_laplace_point_potential( + point_reference, + density_mass, scale_factor, ) remainders = tuple( @@ -395,16 +321,16 @@ def point_density(r: float, t: float) -> float: return NearfieldTemplateExperiment( fan=fan, order=order, - source_order=source_order, near_point_target=near_point_target, point_target_reference=point_reference, point_target_low_order=point_low_order, point_target_abs_error=abs(point_low_order - point_reference), - self_interaction=self_value, + scaled_near_point_target=scaled_target, + scaled_point_target_potential=scaled_value, + expected_scaled_point_target_potential=expected_scaled, + scaled_abs_error=abs(scaled_value - expected_scaled), signed_area=fan.signed_area, + density_mass=density_mass, scale_factor=scale_factor, - scaled_self_interaction=scaled_value, - expected_scaled_self_interaction=expected_scaled, - scaled_abs_error=abs(scaled_value - expected_scaled), diagonal_remainders=remainders, ) diff --git a/tests/evals/test_nearfield_templates.py b/tests/evals/test_nearfield_templates.py index 126a9e2..d3c60e1 100644 --- a/tests/evals/test_nearfield_templates.py +++ b/tests/evals/test_nearfield_templates.py @@ -5,10 +5,9 @@ from cutkit.evals import ( FanTemplateMap2D, diagonal_remainder_sample, - expected_scaled_laplace_self_interaction, + expected_scaled_laplace_point_potential, point_target_laplace_potential, run_nearfield_template_experiment, - self_interaction_laplace, template_density_mass, ) from cutkit.quadrature import gauss_legendre_01 @@ -28,40 +27,30 @@ def test_fan_map_area_and_metric_are_consistent() -> None: assert metric[0][1] == pytest.approx(metric[1][0]) -def test_self_interaction_obeys_log_kernel_scale_law() -> None: +def test_point_target_potential_obeys_log_kernel_scale_law() -> None: fan = FanTemplateMap2D( vertex=(0.0, 0.0), edge_start=(1.0, 0.0), edge_end=(0.35, 0.9), ) scale_factor = 2.25 - value = self_interaction_laplace(fan, order=5, source_order=6) - scaled = self_interaction_laplace( + target = (0.8, 0.7) + value = point_target_laplace_potential(fan, target, order=6) + scaled = point_target_laplace_potential( fan.scaled(scale_factor), - order=5, - source_order=6, + (scale_factor * target[0], scale_factor * target[1]), + order=6, ) - expected = expected_scaled_laplace_self_interaction( + expected = expected_scaled_laplace_point_potential( value, - fan.signed_area, + template_density_mass(fan, order=6), scale_factor, ) assert scaled == pytest.approx(expected, abs=1.0e-13) -def test_self_interaction_rejects_overlapping_template_nodes() -> None: - fan = FanTemplateMap2D( - vertex=(0.0, 0.0), - edge_start=(1.0, 0.0), - edge_end=(0.35, 0.9), - ) - - with pytest.raises(ValueError, match="must not overlap"): - self_interaction_laplace(fan, order=3, source_order=5) - - -def test_scaled_self_interaction_uses_density_weighted_masses() -> None: +def test_scaled_point_target_potential_uses_density_weighted_mass() -> None: fan = FanTemplateMap2D( vertex=(0.0, 0.0), edge_start=(1.0, 0.0), @@ -71,29 +60,24 @@ def test_scaled_self_interaction_uses_density_weighted_masses() -> None: def source_density(r: float, t: float) -> float: return 1.0 + r - def target_density(r: float, t: float) -> float: - return 1.0 - 0.25 * t - scale_factor = 1.4 - value = self_interaction_laplace( + target = (0.8, 0.7) + value = point_target_laplace_potential( fan, + target, order=4, - source_order=5, source_density=source_density, - target_density=target_density, ) - scaled = self_interaction_laplace( + scaled = point_target_laplace_potential( fan.scaled(scale_factor), + (scale_factor * target[0], scale_factor * target[1]), order=4, - source_order=5, source_density=source_density, - target_density=target_density, ) - expected = expected_scaled_laplace_self_interaction( + expected = expected_scaled_laplace_point_potential( value, template_density_mass(fan, order=4, density=source_density), scale_factor, - target_mass=template_density_mass(fan, order=4, density=target_density), ) assert scaled == pytest.approx(expected, abs=1.0e-13) @@ -163,6 +147,10 @@ def test_baseline_experiment_records_feasibility_checks() -> None: assert result.scaled_abs_error < 1.0e-13 assert result.point_target_abs_error > 0.0 + assert result.scaled_point_target_potential == pytest.approx( + result.expected_scaled_point_target_potential, + abs=1.0e-13, + ) assert result.diagonal_remainders[-1].delta == pytest.approx(1.0e-3) assert abs(result.diagonal_remainders[-1].remainder) < abs( result.diagonal_remainders[0].remainder From f4363a270fbc1a8d380dba6168579b4a233bac4c Mon Sep 17 00:00:00 2001 From: xywei Date: Fri, 12 Jun 2026 19:56:46 -0500 Subject: [PATCH 07/15] docs: use high-order point-target singular split Revise the point-target singular split to use a high-order Taylor-jet distance model instead of a leading metric-only model. Also adjust display equations to use aligned environments and scalable matrix delimiters so GitHub Markdown renders the formulas correctly.\n\nValidation: uv run python scripts/check_docs_freshness.py Entire-Checkpoint: d6bbf8d080ee --- docs/nearfield-template-experiments.md | 176 ++++++++++++++++--------- 1 file changed, 116 insertions(+), 60 deletions(-) diff --git a/docs/nearfield-template-experiments.md b/docs/nearfield-template-experiments.md index 5ef90e3..3439f9b 100644 --- a/docs/nearfield-template-experiments.md +++ b/docs/nearfield-template-experiments.md @@ -11,16 +11,20 @@ folded decomposition: a fan from a seed point `V` to a smooth Bezier trim curve `C(t)`. A degree-`p` Bezier edge is $$ -C(t) = \sum_{a=0}^p B_a^p(t)P_a, -\qquad B_a^p(t)=\binom{p}{a}(1-t)^{p-a}t^a, +\begin{aligned} +C(t) &= \sum_{a=0}^p B_a^p(t)P_a, \\ +B_a^p(t) &= \binom{p}{a}(1-t)^{p-a}t^a, \qquad 0\le t\le 1. +\end{aligned} $$ The folded chart is $$ -T(r,t) = (1-r)V + rC(t), +\begin{aligned} +T(r,t) &= (1-r)V + rC(t), \qquad 0\le r,t\le 1. +\end{aligned} $$ Its Jacobian is @@ -80,72 +84,106 @@ The experiment should classify target/source geometry as: ## Point-Target Singular Split -For a target point near the source chart, let `z_x` be a closest or projected -template coordinate with `T(z_x)` near `x`, and write `xi = z_x + delta`. Then +For a target point near the source chart, choose a template coordinate `z_x` +whose mapped source point `T(z_x)` is closest to `x`, or otherwise represents the +local source point responsible for the near-singular behavior. Write nearby +source coordinates as $$ -T(z_x+\delta)-T(z_x) = DT(z_x)\delta + O(|\delta|^2), +\begin{aligned} +\xi &= z_x + \delta, \\ +d &= x - T(z_x). +\end{aligned} $$ +Here `delta` is the local source displacement away from the nearest source point, +and `d` is the physical target offset from that point. Instead of keeping only +the first metric term, use a high-order Taylor jet of the chart around `z_x`: + $$ -|T(z_x+\delta)-T(z_x)| = -\sqrt{\delta^T M(z_x)\delta}\,(1+O(|\delta|)), -\qquad M(z_x) = DT(z_x)^TDT(z_x). +\begin{aligned} +T(z_x+\delta) +&= T(z_x) + \sum_{1\le |\alpha|\le K} \\ +&\quad \frac{1}{\alpha!}\partial^\alpha T(z_x)\delta^\alpha \\ +&\quad + R_{K+1}(\delta). +\end{aligned} $$ -For an on-surface or asymptotically close target, this gives the leading model +Define the order-`K` displacement polynomial $$ -G(x,T(z_x+\delta)) -= -\frac{1}{2\pi}\log\sqrt{\delta^T M(z_x)\delta} -+ \text{lower-order correction}. +\begin{aligned} +P_K(\delta;z_x,d) +&= d - \sum_{1\le |\alpha|\le K} \\ +&\quad \frac{1}{\alpha!}\partial^\alpha T(z_x)\delta^\alpha. +\end{aligned} $$ -For an off-surface target, include the normal/offset residual `d = x-T(z_x)`: +The singular model for the point-target kernel is then $$ -|x-T(z_x+\delta)|^2 -= |d-DT(z_x)\delta|^2 + \text{higher-order chart terms}. +\begin{aligned} +G_K^{\mathrm{sing}}(x,z_x,\delta) +&= -\frac{1}{2\pi}\log |P_K(\delta;z_x,d)|. +\end{aligned} $$ -The reusable part would be template-space singular model integrals or correction -operators. The runtime payload remains map coefficients, Jacobian data, source -coefficients, target-node offsets, and correction moment/interpolation data. For -a Bezier fan map, `T(r,t)` is polynomial in `(r,t)` and contains mixed -higher-order terms from `rC(t)`. Subtracting only the local metric model leaves -an `O(|delta|)` correction whose first derivative at the singular point can -depend on approach direction. A smooth-remainder expansion should therefore -either include higher-order distance terms in the singular model or treat this -first prototype as a leading singular split plus a bounded local correction. +For an on-surface target, `d=0`. For an off-surface target, `d` carries the +normal and tangential target offset. Increasing `K` moves curvature and mixed +terms from the remainder into the singular model. If `C(t)` is a degree-`p` +Bezier curve, the fan map `T(r,t)=(1-r)V+rC(t)` is a polynomial of total degree +`p+1` in `(r,t)`, so the Taylor jet terminates once `K >= p+1`. For +rational/NURBS-derived charts, `K` is a truncation order chosen by the requested +accuracy. + +The reusable part is a family of template-space singular model integrals or +correction operators parameterized by the finite jet data +$\{\partial^\alpha T(z_x)\}$ and the target offset `d`. The runtime payload remains +map coefficients, Jacobian data, source coefficients, target-node offsets, and +high-order correction or interpolation data. ## Metric-Field Expansion Tables The table-building objective is to separate fixed singular template integrals from per-chart smooth geometry and per-target offset coefficients. For a target node represented by offset coordinates `tau`, use local source coordinates near -the closest chart point and write `xi = z_x + delta`. The on-surface singular +the closest chart point and write `xi = z_x + delta`. The high-order singular distance model is $$ -|T(z_x+\delta)-T(z_x)|^2 -= \delta^T M(z_x)\delta + \text{higher-order chart terms}, -\qquad M(z_x)=DT(z_x)^TDT(z_x). +|x-T(z_x+\delta)|^2 \approx |P_K(\delta;z_x,d)|^2. +$$ + +For `K=1` and an on-surface target, this reduces to the metric model + +$$ +\begin{aligned} +|P_1(\delta;z_x,0)|^2 &= \delta^T M(z_x)\delta, \\ +M(z_x) &= DT(z_x)^TDT(z_x). +\end{aligned} $$ -Choose a positive-definite reference metric `M0` for one metric bin, chart -family, or local average, and write +The higher-order table strategy uses the full finite jet, not just `M(z_x)`. The +metric field remains the first term and a useful organizing parameter, but the +practical expansion variables are the coefficients of `P_K` plus the target +offset `d`. + +Choose a reference jet, beginning with a positive-definite reference metric `M0` +and reference higher-order coefficients for one metric/curvature bin, chart +family, or local average. For the first metric term, write $$ M(z) = M_0 + \Delta M(z). $$ -Then the logarithmic singular factor can be expanded as +Then the metric contribution to the logarithmic singular factor can be expanded as $$ +\begin{aligned} \log(\delta^T M(z)\delta) -= \log(\delta^T M_0\delta) -+ \log\left(1+ -\frac{\delta^T\Delta M(z)\delta}{\delta^T M_0\delta}\right). +&= \log(\delta^T M_0\delta) \\ +&\quad + \log\left(1+\frac{\delta^T\Delta M(z)\delta}{\delta^T M_0\delta}\right). +\end{aligned} $$ If the metric family is binned or normalized so that @@ -157,34 +195,41 @@ $$ then $$ +\begin{aligned} \log(\delta^T M(z)\delta) -= \log(\delta^T M_0\delta) -+ \sum_{k\ge 1}\frac{(-1)^{k+1}}{k} -\left(\frac{\delta^T\Delta M(z)\delta}{\delta^T M_0\delta}\right)^k. +&= \log(\delta^T M_0\delta) \\ +&\quad + \sum_{k\ge 1}\frac{(-1)^{k+1}}{k}\left(\frac{\delta^T\Delta M(z)\delta}{\delta^T M_0\delta}\right)^k. +\end{aligned} $$ -Thus the kernel singular part has the schematic expansion +The same expansion idea applies to the full high-order polynomial distance by +expanding the coefficients of `P_K` around the reference jet. Thus the kernel +singular part has the schematic expansion $$ -G_{\mathrm{sing}}(x,T(z_x+\delta)) -\approx \sum_\alpha c_\alpha(z_x,\tau)S_\alpha(\delta;M_0), +\begin{aligned} +G_K^{\mathrm{sing}}(x,z_x,\delta) +&\approx \sum_\alpha c_\alpha(z_x,\tau)S_\alpha(\delta;\mathcal J_0). +\end{aligned} $$ -where `S_alpha` are fixed singular template functions for the selected reference -metric and `c_alpha(z_x,tau)` are smooth functions of the entries of -`Delta M(z_x)` and the target offset data. The expansion can be truncated either -by polynomial order in `Delta M`, by interpolation in metric/offset space, or by -a learned/empirical low-rank basis. -For full smooth-remainder tables near a chart diagonal, the singular model should -also include enough higher-order chart-distance terms to remove direction- -dependent local corrections. +where $\mathcal J_0$ denotes the selected reference jet. The fixed functions +`S_alpha` depend only on template displacement and the reference jet, while +`c_alpha(z_x,tau)` are smooth functions of metric entries, higher-order chart +derivatives, and target offset data. The expansion can be truncated by polynomial +order in the jet perturbation, by interpolation in jet/offset space, or by a +learned/empirical low-rank basis. For one target node and one source density expansion mode `rho_j`, the singular contribution becomes $$ -u_j^{\mathrm{sing}}(\tau) \approx \sum_\alpha \int_{[0,1]^2} -\rho_j(\xi)c_\alpha(z_x,\tau)S_\alpha(\xi-z_x;M_0)J(\xi)\,d\xi. +\begin{aligned} +u_j^{\mathrm{sing}}(\tau) +&\approx \sum_\alpha \int_{[0,1]^2} \\ +&\quad \rho_j(\xi)c_\alpha(z_x,\tau) \\ +&\quad \times S_\alpha(\xi-z_x;\mathcal J_0)J(\xi)\,d\xi. +\end{aligned} $$ Expand the smooth per-chart/per-target factor in a template basis `p_beta`: @@ -196,13 +241,18 @@ $$ Then $$ -u_j^{\mathrm{sing}}(\tau) \approx -\sum_{\alpha,\beta}q_{\alpha\beta}(\tau)T_{j\alpha\beta}, +\begin{aligned} +u_j^{\mathrm{sing}}(\tau) +&\approx \sum_{\alpha,\beta}q_{\alpha\beta}(\tau)T_{j\alpha\beta}. +\end{aligned} $$ $$ -T_{j\alpha\beta} = \int_{[0,1]^2} -\rho_j(\xi)p_\beta(\xi)S_\alpha(\xi-z_x;M_0)\,d\xi. +\begin{aligned} +T_{j\alpha\beta} +&= \int_{[0,1]^2} \\ +&\quad \rho_j(\xi)p_\beta(\xi)S_\alpha(\xi-z_x;\mathcal J_0)\,d\xi. +\end{aligned} $$ The tensors `T_{j alpha beta}` are precomputed on fixed source template domains, @@ -220,16 +270,20 @@ $$ the metric entries are explicit: $$ -T_r = C(t)-V, -\qquad T_t = rC'(t), +\begin{aligned} +T_r &= C(t)-V, \\ +T_t &= rC'(t). +\end{aligned} $$ $$ M(r,t) = -\begin{bmatrix} +\left[ +\begin{matrix} |C(t)-V|^2 & r(C(t)-V)\cdot C'(t) \\ r(C(t)-V)\cdot C'(t) & r^2|C'(t)|^2 -\end{bmatrix}. +\end{matrix} +\right]. $$ For Bezier or piecewise-Bezier CAD trims, these entries are smooth @@ -301,8 +355,10 @@ kernel. If both the source fan and physical target point are scaled by `lambda`, then $$ +\begin{aligned} u_\lambda(\lambda x) -= \lambda^2\left(u(x) - \frac{\log(\lambda)m_\rho}{2\pi}\right), +&= \lambda^2\left(u(x) - \frac{\log(\lambda)m_\rho}{2\pi}\right). +\end{aligned} $$ where `m_rho` is the density-weighted signed source mass on the unscaled fan From 5f3f7d69e586547b210dc0fe413c4602af258846 Mon Sep 17 00:00:00 2001 From: xywei Date: Fri, 12 Jun 2026 21:04:00 -0500 Subject: [PATCH 08/15] docs: motivate point-target displacement polynomial Clarify that P_K is the high-order local model of the physical displacement x - T(xi), obtained by substituting the chart Taylor jet into the point-target kernel distance. This makes the critical table object explicit before defining it.\n\nValidation: uv run python scripts/check_docs_freshness.py Entire-Checkpoint: b61a31d972cf --- docs/nearfield-template-experiments.md | 21 ++++++++++++++++++--- 1 file changed, 18 insertions(+), 3 deletions(-) diff --git a/docs/nearfield-template-experiments.md b/docs/nearfield-template-experiments.md index 3439f9b..74a40b1 100644 --- a/docs/nearfield-template-experiments.md +++ b/docs/nearfield-template-experiments.md @@ -97,8 +97,17 @@ d &= x - T(z_x). $$ Here `delta` is the local source displacement away from the nearest source point, -and `d` is the physical target offset from that point. Instead of keeping only -the first metric term, use a high-order Taylor jet of the chart around `z_x`: +and `d` is the physical target offset from that point. The kernel depends on the +physical displacement from the source point to the target: + +$$ +x - T(\xi) = x - T(z_x+\delta). +$$ + +The purpose of the singular split is to approximate this displacement accurately +near `delta=0`, while keeping the approximation expressed on the fixed source +template. Instead of keeping only the first metric term, replace the chart by its +order-`K` Taylor jet around `z_x`: $$ \begin{aligned} @@ -109,7 +118,9 @@ T(z_x+\delta) \end{aligned} $$ -Define the order-`K` displacement polynomial +Substituting this jet into `x - T(z_x+delta)` gives a polynomial approximation to +the target-to-source displacement. This is the critical object to tabulate +against: $$ \begin{aligned} @@ -119,6 +130,10 @@ P_K(\delta;z_x,d) \end{aligned} $$ +Thus `P_K` is not an extra geometric map; it is the high-order local model of +the physical vector `x - T(xi)`. Its coefficients are exactly the chart jet at +`z_x` plus the target offset `d`. + The singular model for the point-target kernel is then $$ From 2908753ba6e076a92fa7fe53924730095ad63942 Mon Sep 17 00:00:00 2001 From: xywei Date: Fri, 12 Jun 2026 21:58:56 -0500 Subject: [PATCH 09/15] docs: clarify full-jet expansion tables Rewrite the metric-field expansion tables section around the full high-order displacement polynomial P_K for both on-surface and off-surface point targets. The metric-only model is now explicitly a K=1, d=0 sanity case rather than the intended table construction.\n\nValidation: uv run python scripts/check_docs_freshness.py Entire-Checkpoint: 26fe5c4309e9 --- docs/nearfield-template-experiments.md | 91 +++++++++++++++----------- 1 file changed, 54 insertions(+), 37 deletions(-) diff --git a/docs/nearfield-template-experiments.md b/docs/nearfield-template-experiments.md index 74a40b1..9ee21ec 100644 --- a/docs/nearfield-template-experiments.md +++ b/docs/nearfield-template-experiments.md @@ -160,16 +160,28 @@ high-order correction or interpolation data. ## Metric-Field Expansion Tables The table-building objective is to separate fixed singular template integrals -from per-chart smooth geometry and per-target offset coefficients. For a target -node represented by offset coordinates `tau`, use local source coordinates near -the closest chart point and write `xi = z_x + delta`. The high-order singular -distance model is +from per-chart smooth geometry and per-target offset coefficients. The table is +not built from the metric term alone. It is built from the high-order +target-to-source displacement polynomial `P_K`. + +For a target node represented by offset coordinates `tau`, write source points +near the closest chart point as `xi = z_x + delta`. The singular distance model is $$ |x-T(z_x+\delta)|^2 \approx |P_K(\delta;z_x,d)|^2. $$ -For `K=1` and an on-surface target, this reduces to the metric model +The coefficients of `P_K` are the runtime geometry payload: + +$$ +\mathcal J_K(z_x,d) += \left(d,\{\partial^\alpha T(z_x):1\le |\alpha|\le K\}\right). +$$ + +This includes both on-surface and off-surface targets. On-surface targets have +`d=0`; off-surface targets have nonzero `d`, which may include both normal and +tangential offset components. The metric field is only the first quadratic part +of this larger jet. In the special case `K=1` and `d=0`, the model reduces to $$ \begin{aligned} @@ -178,48 +190,50 @@ M(z_x) &= DT(z_x)^TDT(z_x). \end{aligned} $$ -The higher-order table strategy uses the full finite jet, not just `M(z_x)`. The -metric field remains the first term and a useful organizing parameter, but the -practical expansion variables are the coefficients of `P_K` plus the target -offset `d`. - -Choose a reference jet, beginning with a positive-definite reference metric `M0` -and reference higher-order coefficients for one metric/curvature bin, chart -family, or local average. For the first metric term, write +That special case is useful for sanity checks, but it is not the intended table +construction. The practical table construction chooses a reference jet +$\mathcal J_0$ for a geometry/target-offset bin and expands the actual jet around +it: $$ -M(z) = M_0 + \Delta M(z). +\mathcal J_K(z_x,d) = \mathcal J_0 + \Delta\mathcal J(z_x,d). $$ -Then the metric contribution to the logarithmic singular factor can be expanded as +Equivalently, write the displacement polynomial as a reference model plus a +smooth perturbation: $$ \begin{aligned} -\log(\delta^T M(z)\delta) -&= \log(\delta^T M_0\delta) \\ -&\quad + \log\left(1+\frac{\delta^T\Delta M(z)\delta}{\delta^T M_0\delta}\right). +P_K(\delta;z_x,d) +&= P_K^0(\delta) + \Delta P_K(\delta;z_x,d). \end{aligned} $$ -If the metric family is binned or normalized so that +Then the log kernel can be expanded around the reference displacement: $$ -\left|\frac{\delta^T\Delta M(z)\delta}{\delta^T M_0\delta}\right| < 1, +\begin{aligned} +\log |P_K(\delta;z_x,d)| +&= \log |P_K^0(\delta)| \\ +&\quad + \frac{1}{2}\log\left(1+R_K(\delta;z_x,d)\right), +\end{aligned} $$ -then +where $$ \begin{aligned} -\log(\delta^T M(z)\delta) -&= \log(\delta^T M_0\delta) \\ -&\quad + \sum_{k\ge 1}\frac{(-1)^{k+1}}{k}\left(\frac{\delta^T\Delta M(z)\delta}{\delta^T M_0\delta}\right)^k. +R_K(\delta;z_x,d) +&= \frac{N_K(\delta;z_x,d)}{|P_K^0(\delta)|^2}, \\ +N_K(\delta;z_x,d) +&= 2P_K^0(\delta)\cdot\Delta P_K(\delta;z_x,d) \\ +&\quad + |\Delta P_K(\delta;z_x,d)|^2. \end{aligned} $$ -The same expansion idea applies to the full high-order polynomial distance by -expanding the coefficients of `P_K` around the reference jet. Thus the kernel -singular part has the schematic expansion +When each bin is chosen so the perturbation ratio is controlled, this expression +can be expanded in powers, interpolation modes, or a low-rank basis in the jet +perturbation coefficients. The resulting singular model has the schematic form $$ \begin{aligned} @@ -228,12 +242,11 @@ G_K^{\mathrm{sing}}(x,z_x,\delta) \end{aligned} $$ -where $\mathcal J_0$ denotes the selected reference jet. The fixed functions -`S_alpha` depend only on template displacement and the reference jet, while -`c_alpha(z_x,tau)` are smooth functions of metric entries, higher-order chart -derivatives, and target offset data. The expansion can be truncated by polynomial -order in the jet perturbation, by interpolation in jet/offset space, or by a -learned/empirical low-rank basis. +Here `S_alpha` are fixed template functions for the selected reference jet, and +`c_alpha(z_x,tau)` are smooth functions of the target offset `d`, the metric +entries, and all higher-order chart derivatives included in `P_K`. For Bezier +fan charts, these jet coefficients are smooth functions of the seed point and +Bezier control points. For one target node and one source density expansion mode `rho_j`, the singular contribution becomes @@ -247,13 +260,17 @@ u_j^{\mathrm{sing}}(\tau) \end{aligned} $$ -Expand the smooth per-chart/per-target factor in a template basis `p_beta`: +If the smooth per-chart/per-target factor is also expanded in a template basis +`p_beta`, $$ -c_\alpha(z_x,\tau)J(\xi) \approx \sum_\beta q_{\alpha\beta}(\tau)p_\beta(\xi). +\begin{aligned} +c_\alpha(z_x,\tau)J(\xi) +&\approx \sum_\beta q_{\alpha\beta}(\tau)p_\beta(\xi), +\end{aligned} $$ -Then +then the runtime evaluation uses precomputed source-template tables: $$ \begin{aligned} @@ -273,7 +290,7 @@ $$ The tensors `T_{j alpha beta}` are precomputed on fixed source template domains, or tabulated over a small target-offset grid. At runtime, each folded chart and target node only supply the coefficients `q_{alpha beta}(tau)` from smooth -metric/Jacobian fields, target offset data, and any higher-order correction or +jet/Jacobian fields, target offset data, and any higher-order correction or remainder representation. For the Bezier fan map From 4f5040c0ba05c14e92a6a2ed42cc433cb01726b6 Mon Sep 17 00:00:00 2001 From: xywei Date: Fri, 12 Jun 2026 22:11:48 -0500 Subject: [PATCH 10/15] docs: update near-field feasibility assessment Revise the feasibility assessment to reflect the full point-target displacement-jet analysis. The reusable object is now framed as reference-jet tables plus compact interpolation or low-rank coefficients, with metric-only organization called out as a possible subcase rather than the main hypothesis.\n\nValidation: uv run python scripts/check_docs_freshness.py Entire-Checkpoint: 9b045011c52a --- docs/nearfield-template-experiments.md | 65 +++++++++++++++----------- 1 file changed, 39 insertions(+), 26 deletions(-) diff --git a/docs/nearfield-template-experiments.md b/docs/nearfield-template-experiments.md index 9ee21ec..6a978ad 100644 --- a/docs/nearfield-template-experiments.md +++ b/docs/nearfield-template-experiments.md @@ -320,30 +320,39 @@ $$ For Bezier or piecewise-Bezier CAD trims, these entries are smooth low-parameter functions of `(r,t)`, the seed `V`, and the Bezier control points -`P_a`. This is the mathematical reason precomputed template tables can apply to -each chart piece: the singular table is fixed after choosing `M0`, while -observed folded-decomposition geometry should enter through a small smooth -expansion of `M`, `J`, and the higher-order correction data. +`P_a`. They are the first terms of the larger displacement-jet payload. The +mathematical reason precomputed template tables may apply is not that `M` alone +is universal, but that the full finite jet $\mathcal J_K$, the Jacobian, and the +target offset vary smoothly across the folded-decomposition chart family. ## Feasibility Result -The answer is qualified yes. Singular or nearly singular folded-piece -interactions can be moved to fixed template domains, but not to a finite exact -table independent of geometry. The singular coordinate form is reusable; the -metric and higher-order geometry terms enter through smooth geometry-dependent -payloads. +The answer is qualified yes. Singular or nearly singular source-folded-piece to +physical-point interactions can be moved to fixed source template domains. The +reusable object is not a finite exact table independent of geometry. It is a +family of reference-jet tables plus interpolation, expansion, or low-rank +coefficients for the runtime displacement jet. -The practical hypothesis is stronger than carrying `M(z)` pointwise at runtime: -because the folded maps have smooth, low-parameter structure away from collapsed -seed faces, the metric field `M(z)` should vary smoothly over the template and -across the geometry families produced by folded decomposition. A functional -expansion in the metric data, for example moments, interpolation coefficients, -or a low-rank basis in `M` and the correction terms, may cover enough practical -cut-piece cases to make reusable near-field corrections worthwhile. +The practical hypothesis is now about the compactness of the full point-target +payload, not only the metric field. For each near target, the relevant runtime +data is -The open numerical question is therefore whether the smooth dependence on `M`, -seed location, and curve coefficients is compact enough to approximate by such a -functional expansion for useful geometry families. +$$ +\mathcal J_K(z_x,d),\quad J(\xi),\quad \rho(\xi),\quad \tau, +$$ + +where $\mathcal J_K$ contains the target offset and all chart derivatives used by +`P_K`. Because Bezier fan maps have smooth, low-parameter structure away from +collapsed seed faces, these jet coefficients should vary smoothly across the +folded-decomposition cases produced by CAD-like trims. A functional expansion in +the jet data, target-offset data, and Jacobian/density factors may therefore +cover enough practical cases to make precomputed near-field tables worthwhile. + +The open numerical question is whether the observed set of jets and target +offsets is compact or low-rank enough after binning by reference jet +$\mathcal J_0$. +If many bins or many expansion modes are required, the technique may not be +worthwhile even though the template formulation is mathematically valid. ## Measurements @@ -354,10 +363,14 @@ For each geometry and interaction case, record: - template singular-correction error; - higher-order correction/remainder approximation error; - dependence on quadrature order; -- dependence on seed location and curve coefficients; -- size, rank, and smoothness of geometry-parameterized correction data; -- how many metric-field expansion modes are needed for the observed folded - decomposition cases. +- dependence on target offset, including on-surface, near-surface, containing-box, + and neighbor-box target nodes; +- dependence on seed location and Bezier curve coefficients; +- size, rank, and smoothness of the reference-jet correction data; +- how many reference-jet bins and jet-expansion modes are needed for the observed + folded-decomposition cases; +- whether metric-only organization is sufficient for any subfamily, or whether + higher-order jet coefficients dominate the correction size. The first geometry family should be CAD-independent but CAD-realistic: quadratic and cubic Bezier fan charts with seed locations chosen to produce positive, @@ -403,8 +416,8 @@ error. The idea is feasible as a CUTKIT experiment, with these limits: - Far-field source clouds can remain ordinary signed quadrature sources. -- Near-field correction reuse should target template singular bases plus - functional expansions in the smooth metric field `M(z)` and higher geometry - terms. +- Near-field correction reuse should target point-target template singular bases + plus functional expansions in the smooth high-order displacement jet, target + offset, Jacobian, and density data. - Volumential should still own tree/list composition; CUTKIT should export the local geometry/operator payloads needed by those lists. From c1264d0f68ab8c28d36c5d1eae1843766cb0e51d Mon Sep 17 00:00:00 2001 From: xywei Date: Fri, 12 Jun 2026 22:19:41 -0500 Subject: [PATCH 11/15] docs: add Gaussian-window near-field split Document the Gaussian-window kernel split proposed for controlling point-target near-field tables. The singular table is now restricted to targets inside or very close to each fan piece, while the smooth remainder is handled by ordinary folded-decomposition quadrature.\n\nValidation: uv run python scripts/check_docs_freshness.py Entire-Checkpoint: 3c186d8f69f9 --- docs/nearfield-template-experiments.md | 87 +++++++++++++++++++++++--- 1 file changed, 77 insertions(+), 10 deletions(-) diff --git a/docs/nearfield-template-experiments.md b/docs/nearfield-template-experiments.md index 6a978ad..93adc26 100644 --- a/docs/nearfield-template-experiments.md +++ b/docs/nearfield-template-experiments.md @@ -293,6 +293,63 @@ target node only supply the coefficients `q_{alpha beta}(tau)` from smooth jet/Jacobian fields, target offset data, and any higher-order correction or remainder representation. +## Gaussian Window Split + +A more controlled near-field experiment is to split the kernel before building +tables. Choose a length scale `sigma`, usually tied to the local box size or fan +diameter, and a smooth radial window `w_sigma(r)` that is near one at `r=0` and +rapidly decays to numerical zero for `r` larger than a few `sigma`. For the 2D +Laplace kernel, write + +$$ +\begin{aligned} +G(x,y) &= G_{\mathrm{sing},\sigma}(x,y) + G_{\mathrm{smooth},\sigma}(x,y), \\ +G_{\mathrm{sing},\sigma}(x,y) &= w_\sigma(|x-y|)G(x,y), \\ +G_{\mathrm{smooth},\sigma}(x,y) &= \left(1-w_\sigma(|x-y|)\right)G(x,y). +\end{aligned} +$$ + +The smooth part can be handled by ordinary folded-decomposition quadrature on +the source fan: + +$$ +\begin{aligned} +u_{\mathrm{smooth},\sigma}(x) +&= \int_{[0,1]^2}G_{\mathrm{smooth},\sigma}(x,T(\xi)) \\ +&\quad \times \rho(T(\xi))J(\xi)\,d\xi. +\end{aligned} +$$ + +The singular-windowed part is numerically local: + +$$ +\begin{aligned} +u_{\mathrm{sing},\sigma}(x) +&= \int_{[0,1]^2}G_{\mathrm{sing},\sigma}(x,T(\xi)) \\ +&\quad \times \rho(T(\xi))J(\xi)\,d\xi. +\end{aligned} +$$ + +Because `G_{sing,sigma}` is compactly supported to numerical tolerance, this +term matters only when the physical target point lies inside the fan piece or +within the chosen window radius of it. For target nodes in neighboring boxes that +are outside this support, the singular table contribution is skipped and the +ordinary folded quadrature of the smooth part is sufficient. + +This changes the table problem substantially. The expensive singular table only +needs to cover a restricted local target set: + +- target points on or inside the source fan; +- target points within a few `sigma` of the fan boundary; +- target offsets whose support intersects the fan under the high-order local + displacement model `P_K`. + +The smooth remainder no longer needs singular quadrature or reference-jet tables; +it is evaluated directly with the same signed folded quadrature machinery used +for far-field source clouds. The open design choices are the window family, the +scale `sigma`, and the criterion used to skip the singular table for targets +outside the window support. + For the Bezier fan map $$ @@ -327,11 +384,13 @@ target offset vary smoothly across the folded-decomposition chart family. ## Feasibility Result -The answer is qualified yes. Singular or nearly singular source-folded-piece to -physical-point interactions can be moved to fixed source template domains. The -reusable object is not a finite exact table independent of geometry. It is a -family of reference-jet tables plus interpolation, expansion, or low-rank -coefficients for the runtime displacement jet. +The answer is more favorable with a Gaussian-window split. Singular or nearly +singular source-folded-piece to physical-point interactions can be moved to fixed +source template domains, and the tabled part can be restricted to targets inside +or very close to the fan piece. The reusable object is still not a finite exact +table independent of geometry. It is a local family of reference-jet tables plus +interpolation, expansion, or low-rank coefficients for the runtime displacement +jet. The practical hypothesis is now about the compactness of the full point-target payload, not only the metric field. For each near target, the relevant runtime @@ -347,12 +406,16 @@ collapsed seed faces, these jet coefficients should vary smoothly across the folded-decomposition cases produced by CAD-like trims. A functional expansion in the jet data, target-offset data, and Jacobian/density factors may therefore cover enough practical cases to make precomputed near-field tables worthwhile. +The Gaussian-window split improves the odds because the table does not need to +represent weakly near or well-separated target interactions; those move to the +smooth folded-quadrature path. The open numerical question is whether the observed set of jets and target offsets is compact or low-rank enough after binning by reference jet $\mathcal J_0$. -If many bins or many expansion modes are required, the technique may not be -worthwhile even though the template formulation is mathematically valid. +If many bins or many expansion modes are required even after windowing, the +technique may not be worthwhile even though the template formulation is +mathematically valid. ## Measurements @@ -361,7 +424,9 @@ For each geometry and interaction case, record: - reference value; - ordinary pulled-forward tensor-product quadrature error; - template singular-correction error; +- Gaussian-window smooth-part quadrature error; - higher-order correction/remainder approximation error; +- sensitivity to the window scale `sigma` and support cutoff; - dependence on quadrature order; - dependence on target offset, including on-surface, near-surface, containing-box, and neighbor-box target nodes; @@ -371,6 +436,8 @@ For each geometry and interaction case, record: folded-decomposition cases; - whether metric-only organization is sufficient for any subfamily, or whether higher-order jet coefficients dominate the correction size. +- how often source-box and neighbor-box target nodes actually require the + singular table after the window-support test. The first geometry family should be CAD-independent but CAD-realistic: quadratic and cubic Bezier fan charts with seed locations chosen to produce positive, @@ -416,8 +483,8 @@ error. The idea is feasible as a CUTKIT experiment, with these limits: - Far-field source clouds can remain ordinary signed quadrature sources. -- Near-field correction reuse should target point-target template singular bases - plus functional expansions in the smooth high-order displacement jet, target - offset, Jacobian, and density data. +- Near-field correction reuse should target a Gaussian-windowed point-target + singular table for targets inside or very close to each fan piece. The smooth + remainder should use ordinary folded-decomposition quadrature. - Volumential should still own tree/list composition; CUTKIT should export the local geometry/operator payloads needed by those lists. From 92d2324e7d054e1adf01aa3d385342f4a27debb1 Mon Sep 17 00:00:00 2001 From: xywei Date: Fri, 12 Jun 2026 22:22:34 -0500 Subject: [PATCH 12/15] docs: add target-centered Duffy near-field path Document the simplification for interior point targets: choose the folded-decomposition anchor at the target and refold the source panel so the logarithmic singularity sits at the Duffy apex. This reduces the table problem to near-boundary/outside targets after Gaussian windowing.\n\nValidation: uv run python scripts/check_docs_freshness.py Entire-Checkpoint: bd206c686530 --- docs/nearfield-template-experiments.md | 45 +++++++++++++++++++++++--- 1 file changed, 41 insertions(+), 4 deletions(-) diff --git a/docs/nearfield-template-experiments.md b/docs/nearfield-template-experiments.md index 93adc26..397745e 100644 --- a/docs/nearfield-template-experiments.md +++ b/docs/nearfield-template-experiments.md @@ -344,6 +344,35 @@ needs to cover a restricted local target set: - target offsets whose support intersects the fan under the high-order local displacement model `P_K`. +There is an additional simplification for targets inside the source cut region. +Folded decomposition does not require a fixed seed; for a point target `x` inside +the cut region, choose the decomposition anchor to be `x` itself and refold the +local source panel around that target. Each resulting fan has the form + +$$ +\begin{aligned} +T_x(r,t) &= (1-r)x + rC(t). +\end{aligned} +$$ + +The point singularity is then exactly at the Duffy apex `r=0`, and the Jacobian +contributes the usual factor `r`. For the logarithmic kernel, the local integrand +has the form + +$$ +\begin{aligned} +G(x,T_x(r,t))J_x(r,t) +&\sim -\frac{1}{2\pi}\log(r|C(t)-x|)\,r\,\det(C(t)-x,C'(t)), +\end{aligned} +$$ + +which is integrable in the Duffy coordinate. This means inside-target cases may +not need a general reference-jet table at all: they can use target-centered +folded/Duffy quadrature as the singular treatment. The remaining table problem is +then concentrated on targets just outside the fan, near boundaries, or otherwise +too close for the smooth quadrature but not eligible for target-centered +refolding. + The smooth remainder no longer needs singular quadrature or reference-jet tables; it is evaluated directly with the same signed folded quadrature machinery used for far-field source clouds. The open design choices are the window family, the @@ -387,10 +416,13 @@ target offset vary smoothly across the folded-decomposition chart family. The answer is more favorable with a Gaussian-window split. Singular or nearly singular source-folded-piece to physical-point interactions can be moved to fixed source template domains, and the tabled part can be restricted to targets inside -or very close to the fan piece. The reusable object is still not a finite exact +or very close to the fan piece. For targets inside the cut region, a +target-centered Duffy refolding may handle the singularity directly, reducing the +need for tabulation even further. The reusable object is still not a finite exact table independent of geometry. It is a local family of reference-jet tables plus interpolation, expansion, or low-rank coefficients for the runtime displacement -jet. +jet, focused on the near-boundary/outside cases that remain after the Duffy and +window-support tests. The practical hypothesis is now about the compactness of the full point-target payload, not only the metric field. For each near target, the relevant runtime @@ -408,7 +440,9 @@ the jet data, target-offset data, and Jacobian/density factors may therefore cover enough practical cases to make precomputed near-field tables worthwhile. The Gaussian-window split improves the odds because the table does not need to represent weakly near or well-separated target interactions; those move to the -smooth folded-quadrature path. +smooth folded-quadrature path. Target-centered refolding improves the odds again +for interior targets because the singularity is placed at the Duffy apex instead +of being represented by a generic local jet table. The open numerical question is whether the observed set of jets and target offsets is compact or low-rank enough after binning by reference jet @@ -438,6 +472,8 @@ For each geometry and interaction case, record: higher-order jet coefficients dominate the correction size. - how often source-box and neighbor-box target nodes actually require the singular table after the window-support test. +- among supported targets, how many can use target-centered Duffy refolding + instead of a reference-jet table. The first geometry family should be CAD-independent but CAD-realistic: quadratic and cubic Bezier fan charts with seed locations chosen to produce positive, @@ -485,6 +521,7 @@ The idea is feasible as a CUTKIT experiment, with these limits: - Far-field source clouds can remain ordinary signed quadrature sources. - Near-field correction reuse should target a Gaussian-windowed point-target singular table for targets inside or very close to each fan piece. The smooth - remainder should use ordinary folded-decomposition quadrature. + remainder should use ordinary folded-decomposition quadrature. Interior targets + should first try target-centered Duffy refolding before falling back to tables. - Volumential should still own tree/list composition; CUTKIT should export the local geometry/operator payloads needed by those lists. From d62167a8e4d3542039fea7bcd6c480acdbc83692 Mon Sep 17 00:00:00 2001 From: xywei Date: Fri, 12 Jun 2026 22:26:11 -0500 Subject: [PATCH 13/15] docs: clarify windowed singular asymptotics Update the Gaussian-window split to reflect the DMK-style approach: evaluate the compact local singular part with analytic or semi-analytic asymptotic expansion, while using ordinary folded quadrature for the smooth remainder. The local table target is now asymptotic moments/coefficient maps rather than a generic near-field table.\n\nValidation: uv run python scripts/check_docs_freshness.py Entire-Checkpoint: 4d7a58573bba --- docs/nearfield-template-experiments.md | 97 ++++++++++++++------------ 1 file changed, 54 insertions(+), 43 deletions(-) diff --git a/docs/nearfield-template-experiments.md b/docs/nearfield-template-experiments.md index 397745e..4799d6e 100644 --- a/docs/nearfield-template-experiments.md +++ b/docs/nearfield-template-experiments.md @@ -293,13 +293,13 @@ target node only supply the coefficients `q_{alpha beta}(tau)` from smooth jet/Jacobian fields, target offset data, and any higher-order correction or remainder representation. -## Gaussian Window Split +## Gaussian Window And Asymptotic Local Part -A more controlled near-field experiment is to split the kernel before building -tables. Choose a length scale `sigma`, usually tied to the local box size or fan -diameter, and a smooth radial window `w_sigma(r)` that is near one at `r=0` and -rapidly decays to numerical zero for `r` larger than a few `sigma`. For the 2D -Laplace kernel, write +A more controlled near-field experiment is to split the kernel before deciding +what must be tabled. Choose a length scale `sigma`, usually tied to the local box +size or fan diameter, and a smooth radial window `w_sigma(r)` that is near one at +`r=0` and rapidly decays to numerical zero for `r` larger than a few `sigma`. +For the 2D Laplace kernel, write $$ \begin{aligned} @@ -333,17 +333,25 @@ $$ Because `G_{sing,sigma}` is compactly supported to numerical tolerance, this term matters only when the physical target point lies inside the fan piece or within the chosen window radius of it. For target nodes in neighboring boxes that -are outside this support, the singular table contribution is skipped and the +are outside this support, the singular local contribution is skipped and the ordinary folded quadrature of the smooth part is sufficient. -This changes the table problem substantially. The expensive singular table only -needs to cover a restricted local target set: +This changes the hard problem substantially. The windowed singular part only +needs a special local treatment on a restricted target set: - target points on or inside the source fan; - target points within a few `sigma` of the fan boundary; - target offsets whose support intersects the fan under the high-order local displacement model `P_K`. +Following the DMK-style strategy, this local treatment should not be ordinary +folded quadrature of the singular kernel. It should use an analytic or +semi-analytic asymptotic expansion of the windowed singular integral. In local +coordinates, the expansion is built from `P_K`, the Gaussian-window scale, and +smooth expansions of the density and Jacobian. Precomputation, if used, should +target reusable asymptotic moments/coefficient maps for this local expansion, +not a generic table for the whole near-field integral. + There is an additional simplification for targets inside the source cut region. Folded decomposition does not require a fixed seed; for a point target `x` inside the cut region, choose the decomposition anchor to be `x` itself and refold the @@ -368,16 +376,16 @@ $$ which is integrable in the Duffy coordinate. This means inside-target cases may not need a general reference-jet table at all: they can use target-centered -folded/Duffy quadrature as the singular treatment. The remaining table problem is -then concentrated on targets just outside the fan, near boundaries, or otherwise -too close for the smooth quadrature but not eligible for target-centered -refolding. +folded/Duffy coordinates together with the local asymptotic/windowed singular +treatment. The remaining difficult cases are targets just outside the fan, near +boundaries, or otherwise too close for the smooth quadrature but not eligible for +target-centered refolding. -The smooth remainder no longer needs singular quadrature or reference-jet tables; -it is evaluated directly with the same signed folded quadrature machinery used -for far-field source clouds. The open design choices are the window family, the -scale `sigma`, and the criterion used to skip the singular table for targets -outside the window support. +The smooth remainder no longer needs singular quadrature or local asymptotics; it +is evaluated directly with the same signed folded quadrature machinery used for +far-field source clouds. The open design choices are the window family, the scale +`sigma`, the asymptotic expansion order, and the criterion used to skip the local +singular treatment for targets outside the window support. For the Bezier fan map @@ -415,14 +423,13 @@ target offset vary smoothly across the folded-decomposition chart family. The answer is more favorable with a Gaussian-window split. Singular or nearly singular source-folded-piece to physical-point interactions can be moved to fixed -source template domains, and the tabled part can be restricted to targets inside -or very close to the fan piece. For targets inside the cut region, a -target-centered Duffy refolding may handle the singularity directly, reducing the -need for tabulation even further. The reusable object is still not a finite exact -table independent of geometry. It is a local family of reference-jet tables plus -interpolation, expansion, or low-rank coefficients for the runtime displacement -jet, focused on the near-boundary/outside cases that remain after the Duffy and -window-support tests. +source template domains, and the local singular treatment can be restricted to +targets inside or very close to the fan piece. For targets inside the cut region, +target-centered Duffy refolding can place the singularity at the apex. The +reusable object is therefore not a broad near-field table. It is more likely a +small set of asymptotic expansion formulas, moments, or coefficient maps for the +windowed local singular integral, plus direct folded quadrature for the smooth +remainder. The practical hypothesis is now about the compactness of the full point-target payload, not only the metric field. For each near target, the relevant runtime @@ -438,18 +445,19 @@ collapsed seed faces, these jet coefficients should vary smoothly across the folded-decomposition cases produced by CAD-like trims. A functional expansion in the jet data, target-offset data, and Jacobian/density factors may therefore cover enough practical cases to make precomputed near-field tables worthwhile. -The Gaussian-window split improves the odds because the table does not need to -represent weakly near or well-separated target interactions; those move to the -smooth folded-quadrature path. Target-centered refolding improves the odds again -for interior targets because the singularity is placed at the Duffy apex instead -of being represented by a generic local jet table. +The Gaussian-window split improves the odds because the local asymptotic +treatment does not need to represent weakly near or well-separated target +interactions; those move to the smooth folded-quadrature path. Target-centered +refolding improves the odds again for interior targets because the singularity is +placed at the Duffy apex instead of being represented by a generic local jet +table. The open numerical question is whether the observed set of jets and target offsets is compact or low-rank enough after binning by reference jet $\mathcal J_0$. -If many bins or many expansion modes are required even after windowing, the -technique may not be worthwhile even though the template formulation is -mathematically valid. +If many asymptotic terms, bins, or local coefficient modes are required even +after windowing and target-centered refolding, the technique may not be +worthwhile even though the template formulation is mathematically valid. ## Measurements @@ -459,21 +467,23 @@ For each geometry and interaction case, record: - ordinary pulled-forward tensor-product quadrature error; - template singular-correction error; - Gaussian-window smooth-part quadrature error; -- higher-order correction/remainder approximation error; +- Gaussian-window singular asymptotic expansion error; +- higher-order correction/remainder approximation error after the asymptotic + local part is removed; - sensitivity to the window scale `sigma` and support cutoff; - dependence on quadrature order; - dependence on target offset, including on-surface, near-surface, containing-box, and neighbor-box target nodes; - dependence on seed location and Bezier curve coefficients; -- size, rank, and smoothness of the reference-jet correction data; -- how many reference-jet bins and jet-expansion modes are needed for the observed - folded-decomposition cases; +- size, rank, and smoothness of the local asymptotic coefficient data; +- how many window/asymptotic orders, bins, and jet modes are needed for the + observed folded-decomposition cases; - whether metric-only organization is sufficient for any subfamily, or whether higher-order jet coefficients dominate the correction size. - how often source-box and neighbor-box target nodes actually require the - singular table after the window-support test. + singular local treatment after the window-support test. - among supported targets, how many can use target-centered Duffy refolding - instead of a reference-jet table. + instead of a more general local expansion. The first geometry family should be CAD-independent but CAD-realistic: quadratic and cubic Bezier fan charts with seed locations chosen to produce positive, @@ -520,8 +530,9 @@ The idea is feasible as a CUTKIT experiment, with these limits: - Far-field source clouds can remain ordinary signed quadrature sources. - Near-field correction reuse should target a Gaussian-windowed point-target - singular table for targets inside or very close to each fan piece. The smooth - remainder should use ordinary folded-decomposition quadrature. Interior targets - should first try target-centered Duffy refolding before falling back to tables. + local asymptotic expansion for targets inside or very close to each fan piece. + The smooth remainder should use ordinary folded-decomposition quadrature. + Interior targets should first try target-centered Duffy refolding before falling + back to the more general local expansion. - Volumential should still own tree/list composition; CUTKIT should export the local geometry/operator payloads needed by those lists. From 5133f90e0abd2c8700645e513209861252d3ef09 Mon Sep 17 00:00:00 2001 From: xywei Date: Mon, 15 Jun 2026 19:28:42 -0500 Subject: [PATCH 14/15] docs: record near-field residual models Add the DMK-style smooth/local split report, local residual model catalogue, high-order smooth-boundary moment notes, and the complement-subtraction singular quadrature route.\n\nExtend the near-field prototype and tests with Ewald log split helpers, local residual moments, curved-boundary/Taylor graph checks, and precomputed-style smooth-boundary interpolation experiments.\n\nThe complement-subtraction route records the identity that a cut-box singular moment can be computed as a full-box singular table lookup minus a target-separated smooth folded complement integral, with boundary and endpoint-node caveats.\n\nTested with uv sync --extra dev uv run prek install --overwrite --hook-type pre-commit --hook-type pre-push Overwriting existing hook at `/home/xywei/Work/fmm/cutkit/.git/hooks/pre-commit` prek installed at `/home/xywei/Work/fmm/cutkit/.git/hooks/pre-commit` Overwriting existing hook at `/home/xywei/Work/fmm/cutkit/.git/hooks/pre-push` prek installed at `/home/xywei/Work/fmm/cutkit/.git/hooks/pre-push` uv run ruff format --check . 88 files already formatted uv run ruff check . All checks passed! uv run mypy src tests Success: no issues found in 75 source files uv run python scripts/check_architecture.py Architecture check passed. uv run python scripts/check_docs_freshness.py Docs freshness check passed. uv run pytest -q ........................................................................ [ 18%] ........................................................................ [ 36%] ........................................................................ [ 54%] .................................................................s...... [ 72%] ........................................................................ [ 90%] ...................................... [100%] uv run python scripts/run_cutpanel_eval.py name | area | cut_fraction | area_err | max_moment_err | status --- | --- | --- | --- | --- | --- unit-square | 1.000000 | 1.000000 | 6.661e-16 | 1.110e-15 | PASS square-with-hole | 0.750000 | 0.750000 | 5.551e-16 | 3.886e-16 | PASS rect-with-thin-hole | 1.980000 | 0.990000 | 1.554e-15 | 3.109e-15 | PASS square-with-seam-adjacent-slot | 0.998200 | 0.998200 | 6.661e-16 | 1.998e-15 | PASS square-with-ultra-thin-frame | 0.011964 | 0.011964 | 1.076e-16 | 5.065e-16 | PASS imported-prod-courtyard-double-cutout | 12.040000 | 0.802667 | 1.066e-14 | 9.237e-14 | PASS imported-prod-long-corridor-offset-room | 11.620000 | 0.806944 | 8.882e-15 | 3.411e-13 | PASS imported-prod-staggered-triple-core | 23.980000 | 0.749375 | 1.776e-14 | 1.023e-12 | PASS imported-prod-offset-courtyard-barbell | 18.360000 | 0.690226 | 1.776e-14 | 3.411e-13 | PASS imported-prod-service-spine-notches | 21.780000 | 0.806667 | 1.066e-14 | 1.933e-12 | PASS fuzz-candidate-dual-slim-slots-a | 0.894400 | 0.894400 | 5.551e-16 | 1.443e-15 | PASS fuzz-candidate-near-seam-pair | 0.958900 | 0.958900 | 5.551e-16 | 6.661e-16 | PASS fuzz-candidate-quad-pocket-micro | 2.686200 | 0.839438 | 1.332e-15 | 5.773e-15 | PASS fuzz-candidate-staggered-slot-grid | 3.445000 | 0.768973 | 3.553e-15 | 7.994e-15 | PASS fuzz-candidate-thin-z-corridor | 2.208000 | 0.766667 | 8.882e-16 | 3.997e-15 | PASS fuzz-candidate-triple-pocket-offset | 2.607000 | 0.846429 | 2.220e-15 | 3.553e-15 | PASS Cut-panel evaluation passed.. Entire-Checkpoint: adf2c0e740b2 --- ...0260612-analytic-local-residual-moments.md | 53 + ...0260612-dmk-split-nearfield-experiments.md | 65 + ...612-local-residual-nearfield-experiment.md | 52 + .../20260613-local-model-catalogue.md | 37 + .../20260613-local-model-moment-tests.md | 64 + ...plement-subtraction-singular-quadrature.md | 49 + docs/index.md | 2 + docs/nearfield-dmk-split-report.md | 193 +++ docs/nearfield-local-model-catalogue.md | 746 +++++++++ docs/nearfield-template-experiments.md | 188 ++- scripts/run_nearfield_template_experiment.py | 147 +- src/cutkit/evals/__init__.py | 44 + src/cutkit/evals/nearfield_templates.py | 1431 ++++++++++++++++- tests/evals/test_nearfield_templates.py | 325 ++++ 14 files changed, 3362 insertions(+), 34 deletions(-) create mode 100644 docs/exec-plans/completed/20260612-analytic-local-residual-moments.md create mode 100644 docs/exec-plans/completed/20260612-dmk-split-nearfield-experiments.md create mode 100644 docs/exec-plans/completed/20260612-local-residual-nearfield-experiment.md create mode 100644 docs/exec-plans/completed/20260613-local-model-catalogue.md create mode 100644 docs/exec-plans/completed/20260613-local-model-moment-tests.md create mode 100644 docs/exec-plans/completed/20260615-complement-subtraction-singular-quadrature.md create mode 100644 docs/nearfield-dmk-split-report.md create mode 100644 docs/nearfield-local-model-catalogue.md diff --git a/docs/exec-plans/completed/20260612-analytic-local-residual-moments.md b/docs/exec-plans/completed/20260612-analytic-local-residual-moments.md new file mode 100644 index 0000000..1360e82 --- /dev/null +++ b/docs/exec-plans/completed/20260612-analytic-local-residual-moments.md @@ -0,0 +1,53 @@ +# Analytic Local Residual Moments + +## Objective + +Test the simplest analytic/asymptotic local-residual approximation for the 2D log +Ewald split: full-space Taylor moments of the compact residual. Use it for both +interior windows and trim-intersecting windows first, then measure where the +boundary error appears. + +## Scope + +- Add a full-space analytic local residual estimator to the near-field split + sweep. +- Compare analytic local residual values against high-order folded local + residual references. +- Report interior, boundary-near, exterior-near, and vertex-near target errors. +- Run the comparison on `ipa` and update the report. + +## Non-Goals + +- Implement boundary-aware moments. +- Implement curved CAD fixtures. +- Replace direct local residual references. + +## Acceptance Criteria + +- [x] Report identifies when full-space analytic moments are accurate enough and + when trim boundaries dominate the error. +- [x] `make dev` passes before handoff. + +## Checklist + +- [x] Add analytic local estimator and report fields. +- [x] Run the comparison on `ipa`. +- [x] Update docs with conclusions. +- [x] Validate with `make dev`. + +## Outcomes + +- The leading full-space moment `sigma^2 rho(x)/4` works for interior windows + when the target is well away from the trim. At `sigma=0.08`, interior analytic + local relative errors ranged from `2.6e-6` to `6.4e-4`. +- The same full-space moment fails for trim-intersecting windows. Median + trim-near analytic local relative errors were `0.31`, `0.79`, and `1.41` for + `sigma=0.08`, `0.16`, and `0.32`, respectively. +- Exterior-near targets are especially bad because the full-space moment includes + material outside the source region. +- The next step is a window/boundary classifier: interior windows use analytic + moments; trim-intersecting windows need boundary-aware corrections. + +## Validation + +Ran `make dev` successfully after implementation and report updates. diff --git a/docs/exec-plans/completed/20260612-dmk-split-nearfield-experiments.md b/docs/exec-plans/completed/20260612-dmk-split-nearfield-experiments.md new file mode 100644 index 0000000..cd7ce7e --- /dev/null +++ b/docs/exec-plans/completed/20260612-dmk-split-nearfield-experiments.md @@ -0,0 +1,65 @@ +# DMK Split Near-Field Experiments + +## Objective + +Run a first numerical check of the preferred near-field direction: split the 2D +log kernel into a smooth Ewald/heat-kernel part plus a localized singular +residual, evaluate the smooth part with folded fan quadrature, and measure when +the residual may require boundary-aware handling. + +## Scope + +- Extend the existing straight-fan near-field prototype with a 2D log Ewald split. +- Compare ordinary folded quadrature convergence for the full kernel versus the + smoothed kernel. +- Compare current robust interior-style seeds against a simple node-barycenter + seed on CAD-like fan fixtures. +- Run the numerical sweep on `ipa` and summarize results in + `docs/nearfield-template-experiments.md` and + `docs/nearfield-dmk-split-report.md`. + +## Non-Goals + +- Implement the final boundary-aware correction tables. +- Add a production near-field API. +- Replace existing folded decomposition seed policy. + +## Acceptance Criteria + +- [x] Experiment script emits machine-readable and Markdown summaries. +- [x] Results include interior, boundary-near, and exterior-near target cases. +- [x] Report states whether the smooth split improves quadrature convergence and + whether seed choice materially changes fan quality in the fixtures. +- [x] `make dev` passes before handoff. + +## Checklist + +- [x] Add Ewald-split utilities and experiment dataclasses. +- [x] Add CLI options for the DMK-style sweep and report output. +- [x] Run the sweep on `ipa`. +- [x] Add a concise report to the near-field notes. +- [x] Validate with `make dev`. + +## Outcomes + +- The 2D log Ewald split made the smoothed folded-quadrature path much easier: + at `order=24`, smooth-part errors were at least about `70x` smaller than direct + full-log errors across the tested cases. +- Median improvement was `1.8e3`, `6.6e2`, and `1.4e3` for `sigma` values + `0.08`, `0.16`, and `0.32`, respectively. +- Boundary-near targets also benefited strongly, so the smooth remainder is not + the bottleneck in these fixtures. +- The node-barycenter seed is viable but not uniformly better than the current + robust interior anchor; seed choice should remain a measured quality parameter. +- The next unresolved issue is the localized singular residual when its window + intersects a rigid trim boundary. + +## Validation + +Ran `make dev` successfully after the implementation and report updates. + +## Risks + +- The straight-fan fixtures may understate CAD curvature effects. +- Boundary-aware residual behavior requires a later local-residual experiment with + curved or trim-intersecting fixtures to be conclusive. diff --git a/docs/exec-plans/completed/20260612-local-residual-nearfield-experiment.md b/docs/exec-plans/completed/20260612-local-residual-nearfield-experiment.md new file mode 100644 index 0000000..3cd5c5b --- /dev/null +++ b/docs/exec-plans/completed/20260612-local-residual-nearfield-experiment.md @@ -0,0 +1,52 @@ +# Local Residual Near-Field Experiment + +## Objective + +Measure the compact singular residual from the DMK-style 2D log split directly, +to decide whether ordinary folded quadrature is adequate or whether analytic or +boundary-aware precomputed handling is needed. + +## Scope + +- Extend the near-field split sweep with local-residual folded quadrature values. +- Report full, smooth, local, and reconstructed errors against high-order + references. +- Run the residual-inclusive sweep on `ipa`. +- Update the near-field DMK report with conclusions. + +## Non-Goals + +- Implement analytic residual moments. +- Implement boundary-aware precomputed correction tables. +- Add curved CAD fan fixtures. + +## Acceptance Criteria + +- [x] Residual-inclusive report identifies whether direct local folded quadrature + is the new bottleneck. +- [x] Results include the same target and seed cases as the smooth split sweep. +- [x] `make dev` passes before handoff. + +## Checklist + +- [x] Add local residual measurements to experiment dataclasses and CLI tables. +- [x] Run the sweep on `ipa`. +- [x] Update reports and notes. +- [x] Validate with `make dev`. + +## Outcomes + +- Direct ordinary folded quadrature of the compact residual is not competitive as + the primary method. +- At `order=24`, smooth-part median relative errors are near `1e-7`, while local + residual median relative errors range from `4.2e-4` to `8.9e-3` depending on + `sigma`. +- Reconstructing `smooth + local` with direct local quadrature reproduces the + direct full-log error, confirming that the residual is the remaining hard part. +- The next implementation choice should be a special local residual method: + analytic/asymptotic moments for interior windows or boundary-aware precomputed + folded-fan corrections for trim-intersecting windows. + +## Validation + +Ran `make dev` successfully after the implementation and report updates. diff --git a/docs/exec-plans/completed/20260613-local-model-catalogue.md b/docs/exec-plans/completed/20260613-local-model-catalogue.md new file mode 100644 index 0000000..c035e39 --- /dev/null +++ b/docs/exec-plans/completed/20260613-local-model-catalogue.md @@ -0,0 +1,37 @@ +# Local Model Catalogue + +## Objective + +Record the local analytic/asymptotic model catalogue for compact near-field +residual corrections in 2D and 3D. + +## Scope + +- Catalog one-feature residual windows: interiors, smooth boundaries, corners, + faces, edges, and vertices. +- Derive the common Taylor/modal moment expansion used by these models. +- Record the 2D log residual radial moments used by the current experiments. +- Link the catalogue from the near-field experiment documentation. + +## Acceptance Criteria + +- The catalogue states the routing rule that each residual support should see at + most one geometric singular feature. +- 2D and 3D local models include moment formulas and leading terms where known. +- The main near-field docs point to the catalogue as the next design reference. +- Documentation validation passes. + +## Checklist + +- [x] Add local model catalogue. +- [x] Add common Taylor/modal expansion. +- [x] Add 2D log residual moment derivation. +- [x] Link from docs index and near-field reports. +- [x] Run validation. + +## Outcome + +Added `docs/nearfield-local-model-catalogue.md`, covering the one-feature local +models and their expansion structure. The next implementation step is a +window/feature classifier followed by the 2D straight-boundary half-plane moment +experiment. diff --git a/docs/exec-plans/completed/20260613-local-model-moment-tests.md b/docs/exec-plans/completed/20260613-local-model-moment-tests.md new file mode 100644 index 0000000..77e32a3 --- /dev/null +++ b/docs/exec-plans/completed/20260613-local-model-moment-tests.md @@ -0,0 +1,64 @@ +# Local Model Moment Tests + +## Objective + +Turn the 2D straight-feature local model catalogue into executable checks for the +log Ewald residual moments. + +## Scope + +- Add a radial moment helper for the 2D local log residual. +- Add a sector monomial moment helper for full-plane, half-plane, and wedge + models at the local feature point. +- Test the helpers against known full-plane values and independent polar + quadrature references. +- Add a practical curved cut-cell check using an exact rational quadratic trim + and folded curve quadrature. +- Add an exact target-centered curved-boundary model with ray clipping and finite + radial Ewald moments. +- Add a high-order Taylor graph boundary model and compare it with the exact + circular boundary model. +- Add a precomputation-style interpolation check for high-order smooth-boundary + moments. + +## Acceptance Criteria + +- Full-plane leading and second moments match closed-form values. +- Half-plane and wedge sector moments match direct numerical quadrature within + the reference quadrature tolerance. +- On a curved cut-cell sample, the half-plane model beats the full-plane model + for a target on the smooth trim. +- The exact curved-boundary model agrees under refinement to better than + `1e-10` relative error. +- The Taylor boundary model reaches below `1e-10` relative error at sufficiently + high geometry order for tested local windows. +- The precomputed table interpolation reaches below `1e-10` relative error on + intermediate local-window samples. +- The near-field targeted tests and full validation pass. + +## Checklist + +- [x] Implement radial moment helper. +- [x] Implement sector monomial moment helper. +- [x] Add full-plane moment tests. +- [x] Add half-plane and wedge reference tests. +- [x] Add curved cut-cell practical test. +- [x] Add machine-precision curved-boundary model refinement test. +- [x] Add high-order Taylor boundary model convergence test. +- [x] Add precomputed boundary-table interpolation test. +- [x] Run full validation. + +## Outcome + +Added 2D log-residual moment helpers for straight local models and tests covering +full-plane, half-plane, and wedge sectors. The numerical reference uses a +log-radius polar quadrature to avoid direct endpoint integration of the +logarithmic singularity. Added a curved rational cut-cell experiment showing that +the half-plane model is only a leading diagnostic for a target on a smooth curved +trim; the full-plane moment incorrectly includes outside-domain mass, while +machine precision requires exact or high-order curved boundary ray clipping with +finite radial moments. The high-order Taylor graph model reaches the `1e-10` +relative target on the rational quarter-circle trim at order 8 for `sigma=0.02` +and at order 6 for `sigma=0.005, 0.01`. A first interpolation-table check over +the normalized boundary scale also reaches the `1e-10` target, while direct +online graph-strip quadrature is not accurate enough near the target endpoint. diff --git a/docs/exec-plans/completed/20260615-complement-subtraction-singular-quadrature.md b/docs/exec-plans/completed/20260615-complement-subtraction-singular-quadrature.md new file mode 100644 index 0000000..738e3ac --- /dev/null +++ b/docs/exec-plans/completed/20260615-complement-subtraction-singular-quadrature.md @@ -0,0 +1,49 @@ +# Complement-Subtraction Singular Quadrature + +## Objective + +Record and validate the strategy + +```text +int_{Omega cap B} K(x,y) p(y) dy = int_B K(x,y) p(y) dy - int_{B \ Omega} K(x,y) p(y) dy +``` + +for cut-box near-field moments, using existing full-box singular tables for the +first term and folded decomposition for the smooth complement term. + +## Scope + +- Document the mathematical route and validity conditions. +- Identify when the complement integral is smooth and when it still needs a + singular/near-singular fallback. +- Connect the route to existing full-box moment tables and folded fan pieces. +- Keep this as a design note first; production implementation comes after the + classifier/table API is settled. + +## Acceptance Criteria + +- `docs/nearfield-template-experiments.md` explains complement subtraction as a + candidate cut-box singular quadrature route. +- `docs/nearfield-local-model-catalogue.md` records the smoothness conditions, + endpoint-node caveats, and fallback cases. +- `docs/index.md` remains current if a new standalone document is added. +- `make dev` passes before handoff. + +## Checklist + +- [x] Add complement-subtraction derivation and algorithm sketch. +- [x] Record smooth-complement conditions and failure modes. +- [x] Validate documentation with `make dev`. +- [x] Move this plan to `docs/exec-plans/completed/` with outcomes. + +## Outcomes + +- Added the full-box-minus-complement identity to + `docs/nearfield-template-experiments.md` as a candidate cut-box singular + quadrature route. +- Recorded the key condition: the complement folded integral is smooth only when + the target is separated from `B \ Omega`. +- Added routing rules and the open-quadrature caveat to + `docs/nearfield-local-model-catalogue.md`, including when to fall back to + boundary moments, wedge/cone models, refinement, or target-centered Duffy. +- Validation passed with `make dev`. diff --git a/docs/index.md b/docs/index.md index b6fe266..c3371d8 100644 --- a/docs/index.md +++ b/docs/index.md @@ -24,6 +24,8 @@ Repository-local documentation is the system of record for CUTKIT. - `docs/volumential-handoff.md`: CUTKIT far/near primitive handoff contract for volumential workflows - `docs/meshmode-cut-overlay.md`: meshmode cut-overlay contract and validation semantics - `docs/nearfield-template-experiments.md`: folded fan near-field template derivation and prototype +- `docs/nearfield-dmk-split-report.md`: DMK-style smooth/local split experiment results +- `docs/nearfield-local-model-catalogue.md`: local residual model catalogue and moment expansions - `docs/poisson-benchmarks.md`: Poisson-oriented benchmark usage and profiles - `docs/poisson-galerkin-solver.md`: immersed Poisson Galerkin solve + validation workflow - `docs/formdsl-parity-benchmarks.md`: shared IGA + DG-SEM formdsl parity benchmark usage diff --git a/docs/nearfield-dmk-split-report.md b/docs/nearfield-dmk-split-report.md new file mode 100644 index 0000000..645fd42 --- /dev/null +++ b/docs/nearfield-dmk-split-report.md @@ -0,0 +1,193 @@ +# DMK-Style Near-Field Split Report + +Date: 2026-06-12 + +Remote runner: `ipa` + +Command: + +```bash +PYTHONPATH=src python3 scripts/run_nearfield_template_experiment.py \ + --dmk-split-report \ + --orders 4,6,8,10,14,18,24 \ + --sigmas 0.08,0.16,0.32 \ + --reference-order 128 +``` + +Raw JSON and Markdown artifacts were generated on `ipa` in the remote run +directory and copied to the local temporary workspace for summarization. A second +residual-inclusive run was generated with the same parameters after adding direct +quadrature of `G_{\mathrm{local},\sigma}`. The raw artifacts are not checked in +because the concise tables below capture the relevant results. + +## Setup + +The experiment uses two straight-edge polygon fixtures, each integrated as a +signed sum of folded fans from a shared seed to each boundary edge. It compares: + +- full 2D log-kernel folded quadrature; +- smooth Ewald/heat-kernel complement folded quadrature; +- the implied localized residual, measured as `full_reference - smooth_reference`. + +The split is + +$$ +G(r)=G_{\mathrm{smooth},\sigma}(r)+G_{\mathrm{local},\sigma}(r), +\qquad +G_{\mathrm{local},\sigma}(r)=\frac{1}{4\pi}E_1(r^2/\sigma^2), +$$ + +where `E1` is the exponential integral. The smoothed complement has a removable +limit at `r=0`, so ordinary folded quadrature no longer sees the log singularity. + +Targets cover interior, boundary-near-inside, boundary-near-outside, and +vertex-near-inside cases. Seed modes are the current robust interior anchor and a +simple node barycenter. + +## Summary + +At the highest tested order, `order=24`, the smoothed folded quadrature was +consistently more accurate than direct full-log quadrature. + +sigma | min improvement | median improvement | max improvement | median local/full +--- | --- | --- | --- | --- +0.08 | 9.252e+01 | 1.805e+03 | 3.986e+04 | 1.347e-02 +0.16 | 7.045e+01 | 6.599e+02 | 5.466e+04 | 4.753e-02 +0.32 | 1.179e+02 | 1.359e+03 | 4.176e+04 | 1.530e-01 + +Target class | min improvement | median improvement | max improvement +--- | --- | --- | --- +boundary_near_inside | 9.327e+01 | 1.158e+04 | 5.466e+04 +boundary_near_outside | 7.045e+01 | 8.439e+02 | 1.479e+04 +interior | 1.479e+02 | 1.415e+03 | 4.176e+04 +vertex_near_inside | 9.366e+01 | 5.099e+02 | 1.890e+03 + +Here `improvement = full_abs_error / smooth_abs_error` against the same +high-order reference. The smallest improvement is still about `70x`; median +improvements are hundreds to tens of thousands depending on target class and +window width. + +Relative errors for the same `order=24` samples show the same pattern: + +sigma | full rel err median | smooth rel err median | local rel err median | reconstructed rel err median +--- | ---: | ---: | ---: | ---: +0.08 | 5.691e-05 | 9.370e-08 | 8.892e-03 | 5.691e-05 +0.16 | 5.691e-05 | 7.368e-08 | 1.629e-03 | 5.691e-05 +0.32 | 5.691e-05 | 6.844e-08 | 4.190e-04 | 5.691e-05 + +Target class | full rel err median | smooth rel err median | local rel err median | reconstructed rel err median +--- | ---: | ---: | ---: | ---: +boundary_near_inside | 1.147e-03 | 7.614e-08 | 1.494e-02 | 1.147e-03 +boundary_near_outside | 4.122e-05 | 5.124e-08 | 1.629e-03 | 4.122e-05 +interior | 4.225e-05 | 4.332e-08 | 8.944e-04 | 4.225e-05 +vertex_near_inside | 1.053e-04 | 1.680e-07 | 1.235e-03 | 1.053e-04 + +The residual-inclusive run evaluates + +$$ +\int G_{\mathrm{local},\sigma}(|x-y|)\rho(y)\,dy +$$ + +by the same ordinary folded quadrature. Its error is essentially the full-log +error, while the reconstructed error from `smooth + local` matches the direct +full-log error. This means the current direct local-residual quadrature has not +solved the near-field problem; it simply moves the difficult part into an +isolated compact term. + +## Seed Comparison + +fixture | seed mode | det ratio | edge-distance ratio +--- | --- | --- | --- +convex_quad | interior_anchor | 1.718e+00 | 1.510e+00 +convex_quad | node_barycenter | 2.220e+00 | 1.550e+00 +skew_quad | interior_anchor | 1.425e+00 | 1.474e+00 +skew_quad | node_barycenter | 1.697e+00 | 1.264e+00 + +The node barycenter is viable on these convex fixtures but is not uniformly +better. It improves the edge-distance ratio on `skew_quad` but worsens the +boundary-determinant ratio on both fixtures. This supports treating the seed as +a measured quality parameter rather than replacing the robust interior anchor +yet. + +## Interpretation + +The split validates the preferred direction for the first experiment: + +- Ordinary folded decomposition is effective for the smoothed kernel. +- Direct folded quadrature of the compact local residual is the new bottleneck; + at `order=24`, local residual median relative error ranges from `4.2e-4` to + `8.9e-3` depending on `sigma`, versus smooth-part median relative errors near + `1e-7`. +- The singularity removal, not seed tuning, is the dominant improvement for the + smooth path in these fixtures. +- Wider windows move more of the full interaction into the localized residual: + median `|local/full|` grows from `1.3e-2` at `sigma=0.08` to `1.5e-1` at + `sigma=0.32`. +- Boundary-near targets still benefit strongly in the smooth path, so the next + question is not whether folded quadrature can handle the smooth remainder; it + can. The remaining question is how to evaluate the localized residual with an + analytic, asymptotic, or precomputed boundary-aware method. + +## Analytic Full-Space Local Moment + +A follow-up run tested the simplest analytic local residual estimator for +`G_{\mathrm{local},\sigma}`: replace the clipped local domain by the full plane +and use the leading Taylor moment + +$$ +u_{\mathrm{local},\sigma}(x) +\approx \frac{\sigma^2}{4}\rho(x). +$$ + +For the current polynomial test density, the isotropic Laplacian correction +vanishes, so this is the natural first full-space analytic estimate. It was used +unchanged for interior, boundary-near, exterior-near, and vertex-near targets. + +At `order=24`, the analytic local relative errors were: + +sigma | case set | analytic local rel err min | median | max | analytic reconstructed rel err median +--- | --- | ---: | ---: | ---: | ---: +0.08 | all | 2.594e-06 | 2.610e-01 | 4.969e+00 | 4.023e-03 +0.08 | interior only | 2.594e-06 | 1.869e-04 | 6.371e-04 | 2.342e-06 +0.08 | trim-near only | 1.892e-01 | 3.145e-01 | 4.969e+00 | 4.891e-03 +0.16 | all | 2.446e-04 | 6.148e-01 | 2.378e+00 | 3.038e-02 +0.16 | interior only | 2.446e-04 | 8.227e-04 | 1.466e-03 | 4.383e-05 +0.16 | trim-near only | 4.213e-01 | 7.880e-01 | 2.378e+00 | 3.997e-02 +0.32 | all | 5.745e-02 | 1.060e+00 | 1.623e+00 | 1.532e-01 +0.32 | interior only | 5.745e-02 | 7.254e-02 | 8.766e-02 | 1.414e-02 +0.32 | trim-near only | 6.560e-01 | 1.409e+00 | 1.623e+00 | 1.866e-01 + +By target class, over all tested `sigma` values and seeds: + +Target class | analytic local rel err min | median | max | analytic reconstructed rel err median +--- | ---: | ---: | ---: | ---: +boundary_near_inside | 1.892e-01 | 4.799e-01 | 7.530e-01 | 2.186e-02 +boundary_near_outside | 1.452e+00 | 2.141e+00 | 4.969e+00 | 5.235e-02 +interior | 2.594e-06 | 9.924e-04 | 8.766e-02 | 4.383e-05 +vertex_near_inside | 2.171e-01 | 7.880e-01 | 1.609e+00 | 3.997e-02 + +This confirms the expected split: + +- Full-space analytic moments work for interior windows when `sigma` is small + compared with distance to the physical boundary. +- The same full-space moments are not acceptable for trim-intersecting windows. + They include outside-domain mass and can be worse than direct local quadrature. +- Larger `sigma` improves smooth-part quadrature but worsens the full-space local + approximation unless the target is far from every trim boundary. + +## Next Experiment + +Add a local-residual method that separates two cases: + +- interior support: use analytic Taylor/radial moments when the window is well + separated from the physical boundary; +- boundary-intersecting support: build boundary-aware moments or a small + precomputed folded-fan correction table. + +The residual-inclusive and analytic-moment runs together indicate that direct +local folded quadrature is not competitive as the primary method, while +full-space analytic moments are only valid for interior windows. The next +implementation should therefore add a window/boundary classifier and route +interior windows to analytic moments, with trim-intersecting windows routed to a +boundary-aware correction experiment. The candidate local models and moment +expansions are catalogued in `docs/nearfield-local-model-catalogue.md`. diff --git a/docs/nearfield-local-model-catalogue.md b/docs/nearfield-local-model-catalogue.md new file mode 100644 index 0000000..f667b78 --- /dev/null +++ b/docs/nearfield-local-model-catalogue.md @@ -0,0 +1,746 @@ +# Near-Field Local Model Catalogue + +This note catalogs the local analytic/asymptotic models for compact residual +kernels after a DMK/Ewald-style split. The guiding rule is: + +> choose the residual support so each local correction sees at most one geometric +> singular feature. + +If a residual window sees multiple unrelated trim arcs, corners, faces, or +vertices, it should be shrunk, the source geometry should be subdivided, or the +case should fall back to a precomputed/direct local correction. + +## Common Expansion Form + +Let the compact residual kernel be radial: + +$$ +K_\sigma(z) = K_{\mathrm{local},\sigma}(|z|), +\qquad z = y-x. +$$ + +The local residual contribution is + +$$ +u_{\mathrm{local},\sigma}(x) += \int_{\Omega_x} K_\sigma(y-x)a(y)\,dy, +$$ + +where `a` denotes density times any smooth geometric/Jacobian factor and +`Omega_x` is the source region inside the effective residual support around the +target. + +Choose a local model feature `F` and local coordinates `z`. Expand the smooth +amplitude at the target or closest feature point: + +$$ +a(x+z) += \sum_{|\alpha|\le K} +\frac{1}{\alpha!}\partial^\alpha a(x)z^\alpha ++ R_{K+1}(z). +$$ + +Then + +$$ +u_{\mathrm{local},\sigma}(x) += \sum_{|\alpha|\le K} +\frac{1}{\alpha!}\partial^\alpha a(x) +M_\alpha(F,\sigma) ++ \mathcal R_{K+1}, +$$ + +with feature moments + +$$ +M_\alpha(F,\sigma) += \int_{\Omega_F} K_\sigma(z)z^\alpha\,dz. +$$ + +For curved features, `Omega_F` is itself expanded around a leading flat, wedge, +or cone model. The resulting moments are products or perturbations of radial +moments and angular moments. + +Equivalently, introduce the scaled coordinate + +$$ +z = \sigma q, +\qquad K_\sigma(z)=\sigma^{-s}k(|q|), +$$ + +where `s` is the kernel scaling power. Then + +$$ +M_\alpha(F,\sigma) += \sigma^{d-s+|\alpha|} +\int_{\widehat\Omega_F} k(|q|)q^\alpha\,dq, +$$ + +up to corrections from finite support and curved geometry. For the 2D log +residual used here, `s = 0`, so the leading mass is `O(sigma^2)`. For 3D kernels +whose residuals scale like `sigma^{-1} k(r/sigma)`, the leading mass is also +`O(sigma^2)`; the exact power should be attached to the selected kernel split. + +In a local modal implementation, write the smooth amplitude as + +$$ +a(x+z) \approx \sum_{\beta\in\mathcal B_K} c_\beta \phi_\beta(z), +$$ + +where the modal basis may be monomials, orthogonal polynomials on the local +feature model, or chart-jet modes. The correction is then + +$$ +u_{\mathrm{local},\sigma}(x) +\approx \sum_{\beta\in\mathcal B_K} c_\beta +\mu_\beta(F,\sigma), +\qquad +\mu_\beta(F,\sigma)=\int_{\Omega_F}K_\sigma(z)\phi_\beta(z)\,dz. +$$ + +The Taylor form below is the monomial version of this modal expansion. + +## Feature Catalogue + +| Dimension | Feature count in support | Leading model | Parameters | +| --- | --- | --- | --- | +| 2D | none | plane | target offset, amplitude jet | +| 2D | one smooth trim | half-plane | signed distance, tangent, curvature jet | +| 2D | one corner/junction | wedge | offset, opening angle, edge curvature jets | +| 3D | none | space | target offset, amplitude jet | +| 3D | one smooth face | half-space | signed distance, frame, principal curvature jet | +| 3D | one edge | dihedral wedge | offset, dihedral angle, face/edge jets | +| 3D | one vertex | polyhedral cone | offset, incident-face cone, edge/face jets | + +The catalogue is complete for one geometric singularity because local CAD/cut +features stratify by codimension: interiors, smooth boundary strata, edges or +corners, and vertices. Multi-feature windows are deliberately excluded from one +model and must be split or handled by fallback. + +## 2D Log Residual Moments + +For the 2D log split used in the first experiments, + +$$ +K_\sigma(r) += \frac{1}{4\pi}E_1(r^2/\sigma^2). +$$ + +Define radial moments + +$$ +R_m(\sigma) += \int_0^\infty K_\sigma(r)r^{m+1}\,dr. +$$ + +Using + +$$ +\int_0^\infty u^q E_1(u)\,du = \frac{\Gamma(q+1)}{q+1}, +$$ + +and `u = r^2/sigma^2`, + +$$ +R_m(\sigma) += \frac{\sigma^{m+2}}{8\pi} +\frac{\Gamma\left(\frac{m}{2}+1\right)}{\frac{m}{2}+1}. +$$ + +In particular, + +$$ +2\pi R_0(\sigma) = \frac{\sigma^2}{4}, +$$ + +which gives the leading full-plane estimate used in the experiments: + +$$ +u_{\mathrm{local},\sigma}(x) +\approx \frac{\sigma^2}{4}a(x). +$$ + +For finite support truncation at `r <= c sigma`, replace `R_m` by + +$$ +R_m(c,\sigma) += \int_0^{c\sigma} K_\sigma(r)r^{m+1}\,dr. +$$ + +These truncated moments can be evaluated with recurrence/incomplete-gamma forms +or precomputed as one-dimensional functions of `c`. + +## 2D Local Models + +### Interior Model + +Condition: the residual support lies inside the source region and away from all +trim features. + +Model domain: + +$$ +\Omega_F = \mathbb R^2. +$$ + +The moments are full angular moments: + +$$ +M_{ij}^{\mathrm{int}} += \int_0^{2\pi}\int_0^\infty +K_\sigma(r)(r\cos\theta)^i(r\sin\theta)^j r\,dr\,d\theta. +$$ + +Odd moments vanish by symmetry. The first terms are + +$$ +\begin{aligned} +u_{\mathrm{int}}(x) +&= 2\pi R_0 a(x) ++ \frac{1}{2}\left(\partial_{xx}a+\partial_{yy}a\right) +\pi R_2 + O(\sigma^6) \\ +&= \frac{\sigma^2}{4}a(x) ++ \frac{\sigma^4}{32}\Delta a(x) ++ O(\sigma^6) +\end{aligned} +$$ + +for the 2D log residual. + +### Smooth Boundary Model + +Condition: the residual support intersects one smooth trim segment and no corner. + +Use local tangent-normal coordinates at the closest trim point `p`: + +$$ +z = \tau t + \eta n, +\qquad x-p = d_t\tau + d_n n. +$$ + +For a flat leading model with the source on `eta >= -d_n`, use target-centered +coordinates + +$$ +z = s(\cos\theta,\sin\theta), +$$ + +with angular/radial domain + +$$ +s\sin\theta \ge -d_n. +$$ + +Equivalently, along a ray `theta`, the lower radial limit is + +$$ +s_0(\theta;d_n)= +\begin{cases} +0, & \sin\theta \ge 0, \\ +\frac{-d_n}{\sin\theta}, & \sin\theta < 0, +\end{cases} +$$ + +for an interior target with `d_n >= 0`. Exterior targets reverse which angular +sector contributes. The moments are + +$$ +M_{ij}^{\mathrm{half}}(d_n,\sigma) += \int_{\Theta(d_n)} +\int_{s_0(\theta)}^\infty +K_\sigma(s)(s\cos\theta)^i(s\sin\theta)^j s\,ds\,d\theta. +$$ + +For `d_n = 0`, this reduces to half-plane moments. The leading term is one half +of the full-plane mass: + +$$ +u_{\mathrm{half}}(x) +\approx \frac{\sigma^2}{8}a(x) +$$ + +when the target is exactly on a straight boundary. + +For curved trims, express the boundary as a graph in tangent-normal coordinates: + +$$ +\eta = h(\tau) += \frac{\kappa}{2}\tau^2 ++ \frac{h_3}{6}\tau^3 ++ \cdots. +$$ + +The source condition becomes + +$$ +d_n + s\sin\theta +\ge h(d_t+s\cos\theta). +$$ + +Expand the boundary perturbation around the flat model. This produces curvature +corrections of the schematic form + +$$ +M_\alpha^{\mathrm{bdry}} += M_\alpha^{\mathrm{half}} ++ \kappa C_{\alpha,1}(d/\sigma)\sigma^{|\alpha|+3} ++ h_3 C_{\alpha,2}(d/\sigma)\sigma^{|\alpha|+4} ++ \cdots. +$$ + +The functions `C` are universal coefficient maps in normalized target offset and +can be tabulated. + +### Corner Model + +Condition: the residual support intersects one polygon corner or one trim +junction, and no other geometric singular feature. + +The leading model is a wedge with opening angle `omega`: + +$$ +\Omega_F = \{(r,\theta): r\ge 0, +\theta_0\le \theta\le \theta_1\}, +\qquad \omega = \theta_1-\theta_0. +$$ + +For a target at the corner, moments are separable: + +$$ +M_{ij}^{\mathrm{wedge}} += R_{i+j}(\sigma) +\int_{\theta_0}^{\theta_1} +(\cos\theta)^i(\sin\theta)^j\,d\theta. +$$ + +The leading mass term is + +$$ +u_{\mathrm{wedge}}(x) +\approx \frac{\omega}{2\pi}\frac{\sigma^2}{4}a(x). +$$ + +For a target offset from the corner, rays have angle-dependent entry and exit +limits from the two boundary rays. The same formula applies with radial limits +`s_0(theta), s_1(theta)`. Curved edges meeting at the corner introduce two graph +perturbations and a corner-angle perturbation series. + +## 3D Local Models + +The same taxonomy applies in 3D. Let + +$$ +z = r\omega, +\qquad \omega\in S^2. +$$ + +Define radial moments + +$$ +R_m^{3D}(\sigma) += \int_0^\infty K_\sigma(r)r^{m+2}\,dr. +$$ + +Then feature moments are radial moments times angular moments over subsets of +the sphere, plus curvature corrections. + +### Interior Model + +Condition: support lies inside the source volume and away from faces, edges, and +vertices. + +Model domain: + +$$ +\Omega_F = \mathbb R^3. +$$ + +Moments are + +$$ +M_\alpha^{\mathrm{int},3D} += R_{|\alpha|}^{3D}(\sigma) +\int_{S^2}\omega^\alpha\,d\omega. +$$ + +Odd moments vanish. The first correction is proportional to the Laplacian of the +amplitude. + +### Smooth Face Model + +Condition: support intersects one smooth boundary face and no edge. + +Use local coordinates `(u,v,n)` at the closest face point. The flat leading model +is a half-space: + +$$ +n \ge -d_n. +$$ + +Moments are spherical cap moments: + +$$ +M_\alpha^{\mathrm{face}} += \int_{S^2}\int_{s_0(\omega)}^\infty +K_\sigma(s)(s\omega)^\alpha s^2\,ds\,d\omega. +$$ + +For a target exactly on a flat face, the leading mass term is half of the +full-space mass. Curvature corrections use the face graph + +$$ +n = h(u,v) += \frac{1}{2}(\kappa_1 u^2 + \kappa_2 v^2) ++ \cdots. +$$ + +The coefficient maps depend on normalized offset `d/sigma` and the principal +curvatures scaled by `sigma`. + +### Edge Model + +Condition: support intersects one edge where two faces meet and no vertex. + +The leading model is a dihedral wedge extruded along the edge tangent. In local +cylindrical coordinates around the edge, + +$$ +z = z_e e + r(\cos\theta n_1 + \sin\theta n_2), +\qquad \theta_0\le\theta\le\theta_1. +$$ + +Equivalently, the angular domain on `S^2` is a lune determined by the two face +planes. Moments are + +$$ +M_\alpha^{\mathrm{edge}} += \int_{\Omega_{\mathrm{dihedral}}} +K_\sigma(|z|)z^\alpha\,dz. +$$ + +For a target on a straight edge, the leading mass is the solid-angle fraction of +the full-space mass: + +$$ +u_{\mathrm{edge}}(x) +\approx \frac{\Omega_{\mathrm{edge}}}{4\pi}M_0^{\mathrm{full}}a(x), +$$ + +where `Omega_edge` is the solid angle of the dihedral wedge. Higher terms depend +on angular monomial moments of the spherical lune. Curved faces and curved edge +tangents add perturbation terms from the two surface graphs and the edge curve. + +### Vertex Model + +Condition: support intersects one vertex and no other independent feature. + +The leading model is a polyhedral cone: + +$$ +\Omega_F = \{r\omega: r\ge 0, +\omega\in A\subset S^2\}. +$$ + +Moments are + +$$ +M_\alpha^{\mathrm{vertex}} += R_{|\alpha|}^{3D}(\sigma) +\int_A \omega^\alpha\,d\omega. +$$ + +The leading mass term is again the solid-angle fraction of the full-space mass: + +$$ +u_{\mathrm{vertex}}(x) +\approx \frac{|A|}{4\pi}M_0^{\mathrm{full}}a(x). +$$ + +For curved CAD vertices, the cone receives perturbations from all incident face +charts and edge curves. These are likely table-driven rather than hand-derived in +early experiments. + +## Routing Rules + +Each residual correction should be routed by the local feature set inside the +effective support radius `R_sigma = c sigma`: + +1. No physical feature in support: use interior moments. +2. Exactly one smooth boundary face or trim segment: use smooth-boundary moments. +3. Exactly one corner in 2D or edge in 3D: use wedge or dihedral moments. +4. Exactly one vertex in 3D: use cone moments. +5. More than one feature or ambiguous classification: shrink `sigma`, subdivide + the source region, or use a direct/precomputed local fallback. + +The normalized parameters for reusable tables are: + +- normalized target offset `d/sigma`; +- opening angle or solid-angle data for wedges/cones; +- curvature coefficients scaled by powers of `sigma`; +- amplitude and Jacobian jets; +- chart orientation and sign. + +## Complement-Subtraction Route + +For cut boxes with an ambient tensor-product density representation, a useful +alternative to direct boundary-aware singular moments is to subtract the outside +of-domain complement from an existing full-box singular table. For a box `B`, +physical region `Omega_B = Omega cap B`, source mode `p_alpha`, and target `x`, + +$$ +\int_{\Omega_B} K(x,y)p_\alpha(y)\,dy += \int_B K(x,y)p_\alpha(y)\,dy +- \int_{B\setminus\Omega}K(x,y)p_\alpha(y)\,dy. +$$ + +The first term is the standard interior-box moment. The complement term can be +evaluated by folded decomposition as an ordinary smooth integral if + +$$ +\operatorname{dist}(x, B\setminus\Omega) > 0. +$$ + +This condition holds for targets strictly inside the cut cell with a positive +distance to the cut boundary. It fails, or becomes numerically fragile, when the +target lies on the physical boundary, when the complement closure contains the +target, or when the target-complement distance is small relative to the requested +accuracy and quadrature order. + +Routing for this route is: + +1. If the target is target-separated from the complement, use full-box table minus + folded complement. +2. If the target lies on one smooth boundary feature, use smooth-boundary moments + or a boundary table instead. +3. If the target lies on an edge, corner, or vertex, route to the corresponding + wedge/cone model or refine. +4. If multiple complement pieces approach the target, shrink the box/window, + subdivide the geometry, or use a target-centered fallback. + +The complement folded quadrature should use open quadrature rules on fan or cone +coordinates. Lobatto endpoint nodes are useful for interpolation and element +interfaces, but they are undesirable for this smooth-complement integral because a +fan endpoint, clipped-boundary endpoint, or coning apex may lie on the target or +on a nearly singular geometric feature. Open rules keep all quadrature nodes in +the smooth interior of each complement piece. + +This route is attractive for List 1 matrix generation because the singular part +is geometry independent and can reuse the same full-box tables as uncut boxes. +The geometry-dependent object is the complement moment matrix, which is smooth +under the separation condition and can be built by folded quadrature or compressed +over cut-geometry parameters. + +## Precomputed Smooth-Boundary Moment Scheme + +For a target on or near one smooth boundary feature, introduce local +tangent-normal coordinates at the closest boundary point: + +$$ +y = x + s\tau + n\nu, +$$ + +where `tau` is tangent and `nu` is the inward normal. The flat model is the +half-plane `n >= 0`. The curved boundary is represented as a local graph + +$$ +n = h(s) += c_2 s^2 + c_3 s^3 + c_4 s^4 + \cdots. +$$ + +The local source domain is + +$$ +n \ge h(s). +$$ + +For a compact residual kernel `K_sigma`, the local correction is + +$$ +I_\sigma += \int_{n\ge h(s)} K_\sigma(s,n)a(s,n)\,ds\,dn, +$$ + +where `a` is density times Jacobian and any smooth chart factor. Rewrite this as +a flat half-plane moment minus the curved strip excluded by the boundary: + +$$ +I_\sigma += \int_{n\ge 0} K_\sigma(s,n)a(s,n)\,ds\,dn +- \int_{0\le n\le h(s)}K_\sigma(s,n)a(s,n)\,ds\,dn. +$$ + +This identity is the conceptual split. The final implementation should not do a +naive online tensor quadrature over the strip; the logarithmic endpoint at the +target makes that too slow for machine precision. Instead, normalize by the +residual scale: + +$$ +s = \sigma q, +\qquad n = \sigma p. +$$ + +Then the scaled boundary is + +$$ +p = H(q;\lambda) += \lambda_2 q^2 + \lambda_3 q^3 + \lambda_4 q^4 + \cdots, +$$ + +with normalized geometry coefficients + +$$ +\lambda_m = c_m\sigma^{m-1}. +$$ + +For the 2D log residual, `K_sigma(s,n) = k(q,p)` after scaling, and + +$$ +I_\sigma += \sigma^2 +\int_{p\ge H(q;\lambda)} k(q,p) +a(x+\sigma(q\tau+p\nu))\,dq\,dp. +$$ + +Expand the smooth amplitude in the same local coordinates: + +$$ +a(x+\sigma(q\tau+p\nu)) += \sum_{i+j\le A} a_{ij}\sigma^{i+j}q^i p^j ++ O(\sigma^{A+1}). +$$ + +The correction becomes a finite contraction: + +$$ +I_\sigma +\approx +\sum_{i+j\le A}a_{ij}\sigma^{2+i+j} +M_{ij}(\lambda), +$$ + +where the normalized smooth-boundary moments are + +$$ +M_{ij}(\lambda) += \int_{p\ge H(q;\lambda)} k(q,p)q^i p^j\,dq\,dp. +$$ + +These `M_ij(lambda)` are the reusable objects. They can be represented in either +of two equivalent ways: + +1. tabulate `M_ij(lambda)` over normalized boundary-jet parameters and + interpolate; +2. expand the moment function itself in geometry parameters: + +$$ +M_{ij}(\lambda) += M_{ij}^{\mathrm{half}} ++ \lambda_2 M_{ij}^{(2)} ++ \lambda_3 M_{ij}^{(3)} ++ \lambda_2^2 M_{ij}^{(22)} ++ \cdots. +$$ + +The current prototype tests the first representation on a rational quarter-circle +trim. Seven precomputed samples in the normalized scale interpolate intermediate +smooth-boundary moment values below `1e-10` relative error. + +### Offline Table Generation + +For each kernel split, expansion order, and model family: + +1. choose a normalized parameter box for `lambda`; +2. choose a monomial or modal basis `q^i p^j` up to amplitude order `A`; +3. for each table node in `lambda`, evaluate `M_ij(lambda)` using a robust + high-order method such as exact ray clipping and finite radial moments; +4. store interpolation coefficients or a compressed modal representation; +5. validate the table against direct high-order local references on curved CAD + fixtures. + +The expensive ray-exit solves belong in this offline generation step, not in the +main runtime path. + +### Runtime Evaluation + +For each target/source local residual interaction: + +1. classify the support as interior, one smooth boundary, corner/wedge, or + fallback; +2. if smooth-boundary, find the closest local boundary point and frame + `(tau, nu)`; +3. compute the scaled boundary coefficients `lambda_m = c_m sigma^(m-1)`; +4. compute the amplitude jet coefficients `a_ij` for density, Jacobian, and chart + factors; +5. interpolate or evaluate the precomputed moments `M_ij(lambda)`; +6. return the contraction + +$$ +\sum_{i+j\le A}a_{ij}\sigma^{2+i+j}M_{ij}(\lambda). +$$ + +For higher dimensions, the same idea applies with local coordinates split into +tangential variables and one normal variable. Smooth faces use graph moments; +edges use wedge/dihedral graph moments; vertices use cone-sector moments. + +## Experiment Order + +The recommended implementation sequence is: + +1. 2D interior full-space moments. Done for the leading term. +2. 2D straight-boundary half-plane moments. +3. 2D polygon-corner wedge moments. +4. 2D curved-boundary perturbations. +5. 3D face half-space moments. +6. 3D edge dihedral moments. +7. 3D vertex cone moments. + +At each stage, compare the analytic model against the high-order local residual +reference and only add precomputed tables where closed-form or low-dimensional +moment evaluation is not sufficient. + +## Implemented Checks + +The first executable checks cover straight-feature 2D moments for the log Ewald +residual: + +- full-plane mass and second moments; +- boundary-point half-plane sector moments; +- corner-point wedge sector moments. + +The implementation exposes the radial moment and angular-sector monomial moment +helpers in `cutkit.evals.nearfield_templates`. Tests compare closed-form sector +moments against direct polar quadrature of the local residual kernel. + +A second practical check uses an exact rational quadratic quarter-circle cut cell +with one curved trim and two straight cell faces. For targets on the smooth curved +trim and `sigma = 0.02, 0.04`, the leading flat half-plane local model is only a +zeroth-order diagnostic: it is percent-level accurate, while the full-plane model +is worse than `100%` relative error because it includes outside-domain mass. + +The production-quality boundary-aware model must include the local curved +geometry. For the same rational curved cut cell, the implemented target-centered +model clips each ray against the exact curved cell and evaluates the radial Ewald +moments analytically along each ray. The remaining angular integral is adaptively +resolved, and the test requires refinement agreement below `1e-10` relative +error. This is the practical conclusion: flat half-plane moments are useful for +routing and leading-order estimates, but machine precision requires exact or +high-order curved boundary geometry plus finite radial moments. + +The first high-order boundary-model experiment confirms this route on the same +smooth rational trim. The local circular boundary is replaced by an order-`K` +Taylor graph in tangent-normal coordinates. Along each target-centered ray, the +model solves for the positive graph exit distance and uses finite radial Ewald +moments. For `sigma = 0.02`, the relative error against the exact circular +boundary model decreases from `O(1e-6)` at quadratic order to `O(1e-11)` at order +6 and below `1e-10` at order 8. For smaller windows `sigma = 0.005, 0.01`, order +6 is already below `1e-10` relative error. + +The first precomputation-style check tabulates the high-order smooth-boundary +moment value as a function of the normalized local boundary scale. Seven table +samples over `sigma = 0.004..0.028` interpolate intermediate samples at +`sigma = 0.010, 0.018, 0.026` with maximum relative error below `1e-10`. This is +not yet the final modal table format, but it shows that the online ray solve can +move into an offline/table-generation phase. A direct tensor graph-strip +quadrature was also tested and is not sufficient for machine precision near the +target logarithmic endpoint; the precomputed object should be a normalized moment +or coefficient table, not a naive online strip quadrature. diff --git a/docs/nearfield-template-experiments.md b/docs/nearfield-template-experiments.md index 4799d6e..adc1c3c 100644 --- a/docs/nearfield-template-experiments.md +++ b/docs/nearfield-template-experiments.md @@ -345,17 +345,87 @@ needs a special local treatment on a restricted target set: displacement model `P_K`. Following the DMK-style strategy, this local treatment should not be ordinary -folded quadrature of the singular kernel. It should use an analytic or -semi-analytic asymptotic expansion of the windowed singular integral. In local -coordinates, the expansion is built from `P_K`, the Gaussian-window scale, and +folded quadrature of the singular kernel. Jiang and Greengard's continuous DMK +construction splits the finest-level local interaction into Gaussian-sum +convolutions, accelerated by separation of variables, plus a highly localized +correction handled by Taylor/asymptotic expansion of the source density and +low-dimensional radial or spherical moments. For CUTKIT, the analogous local +coordinates should build those moments from `P_K`, the Gaussian-window scale, and smooth expansions of the density and Jacobian. Precomputation, if used, should target reusable asymptotic moments/coefficient maps for this local expansion, not a generic table for the whole near-field integral. -There is an additional simplification for targets inside the source cut region. -Folded decomposition does not require a fixed seed; for a point target `x` inside -the cut region, choose the decomposition anchor to be `x` itself and refold the -local source panel around that target. Each resulting fan has the form +One important difference from the continuous DMK setting is the role of physical +boundaries. DMK's continuous-source derivation treats smooth densities on a +rectangular box, so artificial leaf-box boundaries do not introduce local jumps +and physical boundary effects are compatible with smooth tapering or extension of +the density. A CUTKIT cut boundary is rigid: the effective density is multiplied +by a domain indicator and can jump across the trim. The first experiment should +therefore try the DMK-inspired split directly in folded-decomposition +coordinates: evaluate the smoothed kernel with ordinary folded quadrature and +evaluate the localized correction only for targets whose window support reaches a +source fan. Measurements should then decide whether a special boundary-aware +local expansion is actually required for windows intersecting the cut boundary, +or whether direct local folded quadrature is adequate at the requested +accuracies. + +### Full-Box Minus Complement Route + +A complementary route is available for cut boxes when the source density is +represented on the ambient box basis. Let `B` be a source leaf box and let +`Omega_B = Omega cap B`. For a target `x in Omega_B` and a local source mode +`p_j(y)`, the desired cut-box moment is + +$$ +I_j^{\Omega_B}(x) += \int_{\Omega_B}K(x,y)p_j(y)\,dy. +$$ + +If a Volumential-style or boxcode-style singular table already provides the +full-box moment + +$$ +I_j^B(x)=\int_B K(x,y)p_j(y)\,dy, +$$ + +then + +$$ +I_j^{\Omega_B}(x) += I_j^B(x)-I_j^{B\setminus\Omega}(x), +\qquad +I_j^{B\setminus\Omega}(x) +=\int_{B\setminus\Omega}K(x,y)p_j(y)\,dy. +$$ + +The first term contains the point singularity and is handled by the existing +interior-box table. The second term is smooth when the target is separated from +the complement `B \setminus Omega`. In that case, folded decomposition can +integrate the complement with ordinary high-order quadrature, using open rules on +fan coordinates so no quadrature node lands on a boundary endpoint or fan apex +that could coincide with the target. + +This identity turns a singular cut-cell integral into: + +1. one reusable singular full-box table lookup; +2. one smooth folded-complement integral over the outside-of-domain portion of the + same box; +3. a subtraction in the same source basis and target normalization. + +It is especially attractive for local List 1 matrices on cut boxes because the +full-box table is geometry independent and the geometry dependence moves to a +smooth complement payload. It is not a replacement for boundary-singular local +models when the target lies on the cut boundary, the complement closure contains +the target, or the target-complement distance is too small for ordinary folded +quadrature to be efficient. Those cases still need boundary-aware moment tables, +target-centered Duffy fallback, smaller boxes, or a residual window that sees only +one feature. + +Target-centered Duffy refolding is useful as a diagnostic or single-target +fallback, but it should not be treated as the main reusable strategy. Folded +decomposition does not require a fixed seed; for one point target `x` inside the +cut region, one can choose the decomposition anchor to be `x` itself and refold +the local source panel around that target. Each resulting fan has the form $$ \begin{aligned} @@ -374,18 +444,23 @@ G(x,T_x(r,t))J_x(r,t) \end{aligned} $$ -which is integrable in the Duffy coordinate. This means inside-target cases may -not need a general reference-jet table at all: they can use target-centered -folded/Duffy coordinates together with the local asymptotic/windowed singular -treatment. The remaining difficult cases are targets just outside the fan, near -boundaries, or otherwise too close for the smooth quadrature but not eligible for -target-centered refolding. +which is integrable in the Duffy coordinate, but not smooth. High-order accuracy +can still require many radial nodes, and the refolding must be rebuilt for each +target point. This makes target-centered Duffy unsuitable as the primary +near-field template mechanism. The main experiment should instead test whether a +shared DMK-inspired local correction, parameterized by target offsets and chart +jets, can amortize over many physical target nodes. Target-centered Duffy remains +a useful reference calculation for validating the local correction on individual +interior targets. The smooth remainder no longer needs singular quadrature or local asymptotics; it is evaluated directly with the same signed folded quadrature machinery used for far-field source clouds. The open design choices are the window family, the scale `sigma`, the asymptotic expansion order, and the criterion used to skip the local -singular treatment for targets outside the window support. +singular treatment for targets outside the window support. This DMK-style split +is the preferred experiment direction: ordinary folded decomposition handles the +smoothed kernel, while only the highly localized residual needs special local +treatment. For the Bezier fan map @@ -419,17 +494,26 @@ mathematical reason precomputed template tables may apply is not that `M` alone is universal, but that the full finite jet $\mathcal J_K$, the Jacobian, and the target offset vary smoothly across the folded-decomposition chart family. +The seed point should also be treated as an experimental degree of freedom. +CUTKIT's current 2D folded quadrature chooses a deterministic strictly interior +anchor for robustness, but near-field quality may improve with other stable +choices. The first alternative should be simple: use the barycenter of the trim +nodes or sampled curve nodes as `V`, then compare it with the current interior +anchor. Later seed-quality criteria can minimize fan aspect ratio, Jacobian +variation, signed-weight cancellation, or the spread of the local jet payloads. + ## Feasibility Result The answer is more favorable with a Gaussian-window split. Singular or nearly singular source-folded-piece to physical-point interactions can be moved to fixed source template domains, and the local singular treatment can be restricted to -targets inside or very close to the fan piece. For targets inside the cut region, -target-centered Duffy refolding can place the singularity at the apex. The -reusable object is therefore not a broad near-field table. It is more likely a -small set of asymptotic expansion formulas, moments, or coefficient maps for the -windowed local singular integral, plus direct folded quadrature for the smooth -remainder. +targets inside or very close to the fan piece. The preferred decomposition is +therefore not a broad near-field table. It is a DMK-inspired split: direct folded +quadrature for the smoothed kernel and a compact local correction for the +singular residual. When the local window lies in a smooth interior neighborhood, +the correction should be analytically derived from Taylor/asymptotic moments. If +the local window intersects a rigid cut boundary, the correction may need +boundary-aware precomputation over folded fan templates. The practical hypothesis is now about the compactness of the full point-target payload, not only the metric field. For each near target, the relevant runtime @@ -448,19 +532,49 @@ cover enough practical cases to make precomputed near-field tables worthwhile. The Gaussian-window split improves the odds because the local asymptotic treatment does not need to represent weakly near or well-separated target interactions; those move to the smooth folded-quadrature path. Target-centered -refolding improves the odds again for interior targets because the singularity is -placed at the Duffy apex instead of being represented by a generic local jet -table. +Duffy refolding should be used only as a validation baseline or fallback because +it is target-specific and leaves a nonsmooth integrable kernel. The open numerical question is whether the observed set of jets and target offsets is compact or low-rank enough after binning by reference jet $\mathcal J_0$. -If many asymptotic terms, bins, or local coefficient modes are required even -after windowing and target-centered refolding, the technique may not be -worthwhile even though the template formulation is mathematically valid. +If many asymptotic terms, bins, boundary-aware tables, or local coefficient modes +are required even after windowing, the technique may not be worthwhile even +though the template formulation is mathematically valid. ## Measurements +The first DMK-style smooth/local split sweep is summarized in +`docs/nearfield-dmk-split-report.md`. It confirms that ordinary folded +decomposition is effective for the smoothed kernel: at `order=24`, the smoothed +folded-quadrature error was at least about `70x` smaller than direct full-log +quadrature error across the tested interior, boundary-near, and vertex-near +targets. The residual-inclusive follow-up shows that direct ordinary folded +quadrature of the compact local residual is the remaining bottleneck: its error +tracks the direct full-log error while the smooth-part error is near `1e-7` +relative. The next unresolved question is therefore which special residual path +to implement first: analytic/asymptotic interior moments or boundary-aware +precomputed folded-fan corrections for trim-intersecting windows. + +The first analytic-moment check used the full-space leading residual moment +`sigma^2 rho(x)/4`. It works for interior windows when the target is well away +from the trim: at `order=24`, interior analytic local relative errors were +`2.6e-6` to `6.4e-4` for `sigma=0.08`. The same formula fails for +trim-intersecting windows, where median analytic local relative errors were +`0.31` for `sigma=0.08`, `0.79` for `sigma=0.16`, and `1.41` for `sigma=0.32`. +The next implementation should therefore classify whether the residual window is +interior to the source region. Interior windows can use analytic moments; +trim-intersecting windows need boundary-aware corrections. The detailed local +model catalogue and moment expansions are recorded in +`docs/nearfield-local-model-catalogue.md`. + +The full-box-minus-complement route should be tested alongside these residual +models. Its expected sweet spot is a cut-box target that is interior to +`Omega cap B` but close enough to the cut boundary that direct folded quadrature +of the singular kernel is poor. The full-box singular table supplies the singular +moment exactly for the ambient basis, and the complement integral should converge +as a smooth folded integral as long as the complement is target-separated. + For each geometry and interaction case, record: - reference value; @@ -474,16 +588,25 @@ For each geometry and interaction case, record: - dependence on quadrature order; - dependence on target offset, including on-surface, near-surface, containing-box, and neighbor-box target nodes; -- dependence on seed location and Bezier curve coefficients; +- dependence on seed location and Bezier curve coefficients, starting with the + current robust interior anchor versus the barycenter of trim or sampled curve + nodes; - size, rank, and smoothness of the local asymptotic coefficient data; - how many window/asymptotic orders, bins, and jet modes are needed for the observed folded-decomposition cases; +- whether the localized singular residual can be handled by analytic interior + moments, or whether boundary-aware precomputed fan-template corrections are + needed when the window intersects a trim boundary; +- whether a full-box singular moment minus a smooth folded complement integral is + more accurate or cheaper than direct boundary-aware residual tables for + target-separated cut-box interactions; - whether metric-only organization is sufficient for any subfamily, or whether higher-order jet coefficients dominate the correction size. - how often source-box and neighbor-box target nodes actually require the singular local treatment after the window-support test. -- among supported targets, how many can use target-centered Duffy refolding - instead of a more general local expansion. +- how much accuracy target-centered Duffy refolding provides as a single-target + baseline for the same node count, and when it is too expensive to use as a + fallback. The first geometry family should be CAD-independent but CAD-realistic: quadratic and cubic Bezier fan charts with seed locations chosen to produce positive, @@ -532,7 +655,10 @@ The idea is feasible as a CUTKIT experiment, with these limits: - Near-field correction reuse should target a Gaussian-windowed point-target local asymptotic expansion for targets inside or very close to each fan piece. The smooth remainder should use ordinary folded-decomposition quadrature. - Interior targets should first try target-centered Duffy refolding before falling - back to the more general local expansion. + The singular residual should use analytic moments where the local window is + interior to the source region; if it intersects a trim boundary, test whether + boundary-aware precomputed folded-fan corrections are necessary. + Target-centered Duffy refolding should be kept as a single-target validation + baseline or fallback, not as the main reusable correction strategy. - Volumential should still own tree/list composition; CUTKIT should export the local geometry/operator payloads needed by those lists. diff --git a/scripts/run_nearfield_template_experiment.py b/scripts/run_nearfield_template_experiment.py index 4e5991f..6fbf473 100755 --- a/scripts/run_nearfield_template_experiment.py +++ b/scripts/run_nearfield_template_experiment.py @@ -4,8 +4,29 @@ from __future__ import annotations import argparse +import json +from dataclasses import asdict +from pathlib import Path -from cutkit.evals import run_nearfield_template_experiment +from cutkit.evals import ( + DmkSplitExperimentReport, + run_dmk_split_nearfield_experiment, + run_nearfield_template_experiment, +) + + +def _parse_int_tuple(value: str) -> tuple[int, ...]: + items = tuple(int(item.strip()) for item in value.split(",") if item.strip()) + if not items: + raise argparse.ArgumentTypeError("expected at least one integer") + return items + + +def _parse_float_tuple(value: str) -> tuple[float, ...]: + items = tuple(float(item.strip()) for item in value.split(",") if item.strip()) + if not items: + raise argparse.ArgumentTypeError("expected at least one float") + return items def _parse_args() -> argparse.Namespace: @@ -22,11 +43,135 @@ def _parse_args() -> argparse.Namespace: default=1.75, help="Positive geometry scale factor for the log-kernel scale-law check.", ) + parser.add_argument( + "--dmk-split-report", + action="store_true", + help="Run the DMK-style smooth/local split experiment instead of the baseline.", + ) + parser.add_argument( + "--orders", + type=_parse_int_tuple, + default=(4, 6, 8, 10, 14, 18), + help="Comma-separated quadrature orders for --dmk-split-report.", + ) + parser.add_argument( + "--sigmas", + type=_parse_float_tuple, + default=(0.08, 0.16, 0.32), + help="Comma-separated Ewald split widths for --dmk-split-report.", + ) + parser.add_argument( + "--reference-order", + type=int, + default=96, + help="Reference quadrature order for --dmk-split-report.", + ) + parser.add_argument( + "--json-output", + type=Path, + help="Optional JSON output path for --dmk-split-report.", + ) + parser.add_argument( + "--markdown-output", + type=Path, + help="Optional Markdown output path for --dmk-split-report.", + ) return parser.parse_args() +def _format_dmk_split_report(report: DmkSplitExperimentReport) -> str: + lines: list[str] = [] + lines.append("# DMK-Style Near-Field Split Report") + lines.append("") + lines.append(f"reference_order: `{report.reference_order}`") + lines.append(f"orders: `{', '.join(str(order) for order in report.orders)}`") + lines.append(f"sigmas: `{', '.join(f'{sigma:g}' for sigma in report.sigmas)}`") + lines.append("") + lines.append("## Seed Quality") + lines.append("") + lines.append( + "fixture | seed_mode | seed | min_det | det_ratio | min_edge_distance | edge_distance_ratio" + ) + lines.append("--- | --- | --- | --- | --- | --- | ---") + for sample in report.seed_quality: + lines.append( + f"{sample.fixture} | {sample.seed_mode} | {sample.seed} | " + f"{sample.min_abs_boundary_det:.6e} | " + f"{sample.max_to_min_abs_boundary_det:.6e} | " + f"{sample.min_edge_distance:.6e} | " + f"{sample.max_to_min_edge_distance:.6e}" + ) + + lines.append("") + lines.append("## Highest-Order Error Summary") + lines.append("") + lines.append( + "fixture | seed_mode | target | sigma | full_error | smooth_error | local_error | analytic_local_error | reconstructed_error | analytic_reconstructed_error | smooth_improvement | local_to_full" + ) + lines.append( + "--- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | ---" + ) + max_order = max(report.orders) + for sample in report.samples: + if sample.order != max_order: + continue + improvement = sample.full_abs_error / max(sample.smooth_abs_error, 1.0e-300) + lines.append( + f"{sample.fixture} | {sample.seed_mode} | {sample.target_label} | " + f"{sample.sigma:.3g} | {sample.full_abs_error:.6e} | " + f"{sample.smooth_abs_error:.6e} | {sample.local_abs_error:.6e} | " + f"{sample.analytic_local_abs_error:.6e} | " + f"{sample.reconstructed_abs_error:.6e} | " + f"{sample.analytic_reconstructed_abs_error:.6e} | " + f"{improvement:.6e} | " + f"{sample.local_to_full_ratio:.6e}" + ) + + lines.append("") + lines.append("## Error By Order") + lines.append("") + lines.append( + "fixture | seed_mode | target | sigma | order | full_error | smooth_error | local_error | analytic_local_error | reconstructed_error | analytic_reconstructed_error | smooth_improvement" + ) + lines.append( + "--- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | ---" + ) + for sample in report.samples: + improvement = sample.full_abs_error / max(sample.smooth_abs_error, 1.0e-300) + lines.append( + f"{sample.fixture} | {sample.seed_mode} | {sample.target_label} | " + f"{sample.sigma:.3g} | {sample.order} | {sample.full_abs_error:.6e} | " + f"{sample.smooth_abs_error:.6e} | {sample.local_abs_error:.6e} | " + f"{sample.analytic_local_abs_error:.6e} | " + f"{sample.reconstructed_abs_error:.6e} | " + f"{sample.analytic_reconstructed_abs_error:.6e} | {improvement:.6e}" + ) + lines.append("") + return "\n".join(lines) + + +def _run_dmk_split_report(args: argparse.Namespace) -> int: + report = run_dmk_split_nearfield_experiment( + orders=args.orders, + sigmas=args.sigmas, + reference_order=args.reference_order, + ) + markdown = _format_dmk_split_report(report) + if args.json_output is not None: + args.json_output.parent.mkdir(parents=True, exist_ok=True) + args.json_output.write_text(json.dumps(asdict(report), indent=2) + "\n") + if args.markdown_output is not None: + args.markdown_output.parent.mkdir(parents=True, exist_ok=True) + args.markdown_output.write_text(markdown) + print(markdown) + return 0 + + def main() -> int: args = _parse_args() + if args.dmk_split_report: + return _run_dmk_split_report(args) + result = run_nearfield_template_experiment( order=args.order, scale_factor=args.scale_factor, diff --git a/src/cutkit/evals/__init__.py b/src/cutkit/evals/__init__.py index 1b9956b..1d3bd92 100644 --- a/src/cutkit/evals/__init__.py +++ b/src/cutkit/evals/__init__.py @@ -76,14 +76,36 @@ ) from cutkit.evals.nearfield_templates import ( DiagonalRemainderSample, + CurvedBoundaryLocalModelReport, + CurvedBoundaryLocalModelSample, + DmkSplitExperimentReport, + DmkSplitSample, FanTemplateMap2D, + GraphBoundaryModelReport, + GraphBoundaryModelSample, NearfieldTemplateExperiment, + PrecomputedBoundaryTableReport, + PrecomputedBoundaryTableSample, + SeedQualitySample, + TaylorBoundaryModelReport, + TaylorBoundaryModelSample, diagonal_remainder_sample, + ewald_log_local_kernel, + ewald_log_local_radial_moment_2d, + ewald_log_local_sector_monomial_moment_2d, + ewald_log_smooth_kernel, expected_scaled_laplace_point_potential, + exponential_integral_e1, laplace_log_kernel, metric_model_distance, point_target_laplace_potential, + point_target_smooth_ewald_potential, + run_curved_boundary_local_model_experiment, + run_dmk_split_nearfield_experiment, + run_graph_boundary_model_experiment, run_nearfield_template_experiment, + run_precomputed_boundary_table_experiment, + run_taylor_boundary_model_experiment, template_density_mass, ) @@ -151,14 +173,36 @@ "run_poisson_galerkin_benchmark", "solve_trimmed_poisson_galerkin", "DiagonalRemainderSample", + "CurvedBoundaryLocalModelReport", + "CurvedBoundaryLocalModelSample", + "DmkSplitExperimentReport", + "DmkSplitSample", "FanTemplateMap2D", + "GraphBoundaryModelReport", + "GraphBoundaryModelSample", "NearfieldTemplateExperiment", + "PrecomputedBoundaryTableReport", + "PrecomputedBoundaryTableSample", + "SeedQualitySample", + "TaylorBoundaryModelReport", + "TaylorBoundaryModelSample", "diagonal_remainder_sample", + "ewald_log_local_kernel", + "ewald_log_local_radial_moment_2d", + "ewald_log_local_sector_monomial_moment_2d", + "ewald_log_smooth_kernel", "expected_scaled_laplace_point_potential", + "exponential_integral_e1", "laplace_log_kernel", "metric_model_distance", "point_target_laplace_potential", + "point_target_smooth_ewald_potential", + "run_curved_boundary_local_model_experiment", + "run_dmk_split_nearfield_experiment", + "run_graph_boundary_model_experiment", "run_nearfield_template_experiment", + "run_precomputed_boundary_table_experiment", + "run_taylor_boundary_model_experiment", "template_density_mass", "FormDslParityBenchmark", "FormDslParityRow", diff --git a/src/cutkit/evals/nearfield_templates.py b/src/cutkit/evals/nearfield_templates.py index 0fdd49f..18e905f 100644 --- a/src/cutkit/evals/nearfield_templates.py +++ b/src/cutkit/evals/nearfield_templates.py @@ -3,14 +3,23 @@ from __future__ import annotations from dataclasses import dataclass -from math import hypot, log, pi, sqrt +from math import atan2, cos, erf, exp, gamma, hypot, isfinite, log, pi, sin, sqrt from typing import Callable -from cutkit.quadrature import gauss_legendre_01 +from cutkit.geometry import ( + CurveEdge2D, + CurveLoop2D, + CurveTrimmedPanel2D, + PanelLoop2D, + TrimmedPanel2D, +) +from cutkit.quadrature import folded_curve_quadrature_rule, gauss_legendre_01 +from cutkit.topology import select_interior_anchor Point2D = tuple[float, float] TemplateDensity2D = Callable[[float, float], float] _COINCIDENT_TOL = 1.0e-14 +_EULER_GAMMA = 0.5772156649015329 def unit_template_density(_r: float, _t: float) -> float: @@ -131,6 +140,161 @@ class NearfieldTemplateExperiment: diagonal_remainders: tuple[DiagonalRemainderSample, ...] +@dataclass(frozen=True) +class DmkSplitSample: + """One quadrature-convergence sample for the DMK-style log-kernel split.""" + + fixture: str + seed_mode: str + seed: Point2D + target_label: str + target: Point2D + sigma: float + order: int + full_value: float + smooth_value: float + local_value: float + analytic_local_value: float + reconstructed_value: float + analytic_reconstructed_value: float + full_reference: float + smooth_reference: float + local_reference: float + full_abs_error: float + smooth_abs_error: float + local_abs_error: float + analytic_local_abs_error: float + reconstructed_abs_error: float + analytic_reconstructed_abs_error: float + local_to_full_ratio: float + + +@dataclass(frozen=True) +class SeedQualitySample: + """Geometry quality summary for one folded-panel seed.""" + + fixture: str + seed_mode: str + seed: Point2D + min_abs_boundary_det: float + max_abs_boundary_det: float + max_to_min_abs_boundary_det: float + min_edge_distance: float + max_edge_distance: float + max_to_min_edge_distance: float + + +@dataclass(frozen=True) +class DmkSplitExperimentReport: + """Report payload for the DMK-style folded near-field split experiment.""" + + reference_order: int + orders: tuple[int, ...] + sigmas: tuple[float, ...] + seed_quality: tuple[SeedQualitySample, ...] + samples: tuple[DmkSplitSample, ...] + + +@dataclass(frozen=True) +class CurvedBoundaryLocalModelSample: + """Practical local-model check on one curved CAD-style cut cell.""" + + sigma: float + order: int + reference_order: int + target: Point2D + folded_value: float + reference_value: float + curved_boundary_value: float + half_plane_value: float + full_plane_value: float + folded_abs_error: float + curved_boundary_abs_error: float + half_plane_abs_error: float + full_plane_abs_error: float + curved_boundary_rel_error: float + half_plane_rel_error: float + full_plane_rel_error: float + + +@dataclass(frozen=True) +class CurvedBoundaryLocalModelReport: + """Report for local residual models on an exact curved cut cell.""" + + sigmas: tuple[float, ...] + orders: tuple[int, ...] + reference_order: int + samples: tuple[CurvedBoundaryLocalModelSample, ...] + + +@dataclass(frozen=True) +class TaylorBoundaryModelSample: + """One high-order Taylor boundary model sample for a smooth trim.""" + + sigma: float + boundary_order: int + exact_value: float + taylor_value: float + abs_error: float + rel_error: float + + +@dataclass(frozen=True) +class TaylorBoundaryModelReport: + """Convergence report for high-order smooth-boundary Taylor models.""" + + sigmas: tuple[float, ...] + boundary_orders: tuple[int, ...] + angular_order: int + samples: tuple[TaylorBoundaryModelSample, ...] + + +@dataclass(frozen=True) +class GraphBoundaryModelSample: + """One graph-strip smooth-boundary model sample.""" + + sigma: float + boundary_order: int + strip_order: int + exact_value: float + graph_value: float + abs_error: float + rel_error: float + + +@dataclass(frozen=True) +class GraphBoundaryModelReport: + """Convergence report for graph-strip boundary corrections.""" + + sigmas: tuple[float, ...] + boundary_orders: tuple[int, ...] + strip_orders: tuple[int, ...] + angular_order: int + samples: tuple[GraphBoundaryModelSample, ...] + + +@dataclass(frozen=True) +class PrecomputedBoundaryTableSample: + """One interpolated precomputed boundary-table sample.""" + + sigma: float + exact_value: float + table_value: float + abs_error: float + rel_error: float + + +@dataclass(frozen=True) +class PrecomputedBoundaryTableReport: + """Report for interpolation over precomputed smooth-boundary moments.""" + + table_sigmas: tuple[float, ...] + eval_sigmas: tuple[float, ...] + boundary_order: int + angular_order: int + samples: tuple[PrecomputedBoundaryTableSample, ...] + + def laplace_log_kernel(target: Point2D, source: Point2D) -> float: """Return the 2D Laplace fundamental solution ``-log(|x-y|)/(2*pi)``.""" @@ -140,6 +304,131 @@ def laplace_log_kernel(target: Point2D, source: Point2D) -> float: return -log(distance) / (2.0 * pi) +def exponential_integral_e1(x: float, *, tol: float = 1.0e-15) -> float: + """Return ``E1(x) = integral_x^inf exp(-t)/t dt`` for ``x > 0``. + + This small stdlib-only evaluator is intended for experiment-scale kernel + splitting, not as a general special-functions package. + """ + + if x <= 0.0 or not isfinite(x): + raise ValueError("x must be positive and finite") + + if x <= 4.0: + term_power = 1.0 + factorial = 1.0 + series = 0.0 + for k in range(1, 256): + term_power *= -x + factorial *= k + term = term_power / (k * factorial) + series += term + if abs(term) <= tol * max(1.0, abs(series)): + break + return -_EULER_GAMMA - log(x) - series + + term = 1.0 + series = 1.0 + previous_abs = abs(term) + for k in range(1, 256): + term *= -k / x + current_abs = abs(term) + if current_abs > previous_abs: + break + series += term + if current_abs <= tol * max(1.0, abs(series)): + break + previous_abs = current_abs + return exp(-x) * series / x + + +def ewald_log_local_kernel(target: Point2D, source: Point2D, *, sigma: float) -> float: + """Return the localized 2D log Ewald residual ``E1(r^2/sigma^2)/(4*pi)``.""" + + if sigma <= 0.0: + raise ValueError("sigma must be positive") + distance = hypot(target[0] - source[0], target[1] - source[1]) + if distance <= 0.0: + raise ValueError("local log kernel is singular at coincident points") + return exponential_integral_e1((distance / sigma) ** 2) / (4.0 * pi) + + +def ewald_log_smooth_kernel(target: Point2D, source: Point2D, *, sigma: float) -> float: + """Return the smooth complement of the localized 2D log Ewald residual.""" + + if sigma <= 0.0: + raise ValueError("sigma must be positive") + distance = hypot(target[0] - source[0], target[1] - source[1]) + if distance <= _COINCIDENT_TOL: + return (_EULER_GAMMA - 2.0 * log(sigma)) / (4.0 * pi) + return laplace_log_kernel(target, source) - ewald_log_local_kernel( + target, + source, + sigma=sigma, + ) + + +def ewald_log_local_radial_moment_2d(*, sigma: float, monomial_degree: int) -> float: + """Return ``int_0^inf K_sigma(r) r^(monomial_degree + 1) dr``. + + Here ``K_sigma(r) = E1(r^2/sigma^2)/(4*pi)`` is the localized 2D log + residual. Multiplying this radial moment by the corresponding angular + monomial moment gives sector moments for full-plane, half-plane, and wedge + local models. + """ + + if sigma <= 0.0: + raise ValueError("sigma must be positive") + if monomial_degree < 0: + raise ValueError("monomial_degree must be nonnegative") + + half_degree = 0.5 * monomial_degree + return ( + sigma ** (monomial_degree + 2) + * gamma(half_degree + 1.0) + / (8.0 * pi * (half_degree + 1.0)) + ) + + +def ewald_log_local_sector_monomial_moment_2d( + *, + sigma: float, + x_power: int, + y_power: int, + theta_start: float, + theta_end: float, + angular_order: int = 96, +) -> float: + """Return a 2D local residual monomial moment over an angular sector. + + The sector has vertex at the target and angle range ``[theta_start, + theta_end]``. This is the leading straight-feature model for full-plane, + half-plane, and wedge residual corrections. + """ + + if x_power < 0 or y_power < 0: + raise ValueError("monomial powers must be nonnegative") + if angular_order < 1: + raise ValueError("angular_order must be positive") + if theta_end < theta_start: + raise ValueError("theta_end must be at least theta_start") + + radial = ewald_log_local_radial_moment_2d( + sigma=sigma, + monomial_degree=x_power + y_power, + ) + if theta_end == theta_start: + return 0.0 + + nodes, weights = gauss_legendre_01(angular_order) + width = theta_end - theta_start + angular = 0.0 + for node, weight in zip(nodes, weights, strict=True): + theta = theta_start + width * node + angular += weight * width * (cos(theta) ** x_power) * (sin(theta) ** y_power) + return radial * angular + + def metric_model_distance( fan: FanTemplateMap2D, *, @@ -218,6 +507,1144 @@ def point_target_laplace_potential( return total +def point_target_smooth_ewald_potential( + fan: FanTemplateMap2D, + target: Point2D, + *, + order: int, + sigma: float, + source_density: TemplateDensity2D | None = None, +) -> float: + """Approximate the smoothed 2D log-kernel source integral at one target.""" + + if order < 1: + raise ValueError("order must be positive") + if source_density is None: + source_density = unit_template_density + + nodes, weights = gauss_legendre_01(order) + total = 0.0 + for r, wr in zip(nodes, weights, strict=True): + jacobian = fan.signed_jacobian(r) + for t, wt in zip(nodes, weights, strict=True): + source = fan.point(r, t) + total += ( + ewald_log_smooth_kernel(target, source, sigma=sigma) + * source_density(r, t) + * jacobian + * wr + * wt + ) + return total + + +def _edge_distance(point: Point2D, start: Point2D, end: Point2D) -> float: + edge = _sub(end, start) + length_sq = _dot(edge, edge) + if length_sq <= 0.0: + return hypot(point[0] - start[0], point[1] - start[1]) + rel = _sub(point, start) + tau = max(0.0, min(1.0, _dot(rel, edge) / length_sq)) + closest = (start[0] + tau * edge[0], start[1] + tau * edge[1]) + return hypot(point[0] - closest[0], point[1] - closest[1]) + + +def _node_barycenter(nodes: tuple[Point2D, ...]) -> Point2D: + scale = 1.0 / len(nodes) + return ( + scale * sum(point[0] for point in nodes), + scale * sum(point[1] for point in nodes), + ) + + +def _panel_fans( + nodes: tuple[Point2D, ...], *, seed: Point2D +) -> tuple[FanTemplateMap2D, ...]: + return tuple( + FanTemplateMap2D( + vertex=seed, + edge_start=nodes[index], + edge_end=nodes[(index + 1) % len(nodes)], + ) + for index in range(len(nodes)) + ) + + +def _panel_seed_quality( + *, + fixture: str, + seed_mode: str, + nodes: tuple[Point2D, ...], + seed: Point2D, +) -> SeedQualitySample: + fans = _panel_fans(nodes, seed=seed) + dets = tuple(abs(fan.boundary_det) for fan in fans) + distances = tuple( + _edge_distance(seed, nodes[index], nodes[(index + 1) % len(nodes)]) + for index in range(len(nodes)) + ) + min_det = min(dets) + min_distance = min(distances) + max_det = max(dets) + max_distance = max(distances) + return SeedQualitySample( + fixture=fixture, + seed_mode=seed_mode, + seed=seed, + min_abs_boundary_det=min_det, + max_abs_boundary_det=max_det, + max_to_min_abs_boundary_det=max_det / min_det + if min_det > 0.0 + else float("inf"), + min_edge_distance=min_distance, + max_edge_distance=max_distance, + max_to_min_edge_distance=( + max_distance / min_distance if min_distance > 0.0 else float("inf") + ), + ) + + +def _physical_density(point: Point2D) -> float: + return 1.0 + 0.25 * point[0] - 0.15 * point[1] + 0.1 * point[0] * point[1] + + +def _rational_quarter_circle_edge() -> CurveEdge2D: + """Return the exact rational quadratic arc used as a curved cut trim.""" + + p0 = (0.0, 0.5) + p1 = (0.0, 0.0) + p2 = (0.5, 0.0) + w0 = 1.0 + w1 = 1.0 / sqrt(2.0) + w2 = 1.0 + + def basis(t: float) -> tuple[float, float, float]: + one_minus_t = 1.0 - t + return (one_minus_t * one_minus_t, 2.0 * one_minus_t * t, t * t) + + def basis_derivative(t: float) -> tuple[float, float, float]: + return (-2.0 * (1.0 - t), 2.0 - 4.0 * t, 2.0 * t) + + def weighted_sum(t: float) -> tuple[float, float, float]: + b0, b1, b2 = basis(t) + denominator = w0 * b0 + w1 * b1 + w2 * b2 + numerator_x = w0 * b0 * p0[0] + w1 * b1 * p1[0] + w2 * b2 * p2[0] + numerator_y = w0 * b0 * p0[1] + w1 * b1 * p1[1] + w2 * b2 * p2[1] + return numerator_x, numerator_y, denominator + + def weighted_sum_derivative(t: float) -> tuple[float, float, float]: + db0, db1, db2 = basis_derivative(t) + denominator_derivative = w0 * db0 + w1 * db1 + w2 * db2 + numerator_x_derivative = w0 * db0 * p0[0] + w1 * db1 * p1[0] + w2 * db2 * p2[0] + numerator_y_derivative = w0 * db0 * p0[1] + w1 * db1 * p1[1] + w2 * db2 * p2[1] + return numerator_x_derivative, numerator_y_derivative, denominator_derivative + + def point(t: float) -> Point2D: + numerator_x, numerator_y, denominator = weighted_sum(t) + return (numerator_x / denominator, numerator_y / denominator) + + def tangent(t: float) -> Point2D: + numerator_x, numerator_y, denominator = weighted_sum(t) + dx_num, dy_num, denominator_derivative = weighted_sum_derivative(t) + denominator_sq = denominator * denominator + return ( + (dx_num * denominator - numerator_x * denominator_derivative) + / denominator_sq, + (dy_num * denominator - numerator_y * denominator_derivative) + / denominator_sq, + ) + + return CurveEdge2D(evaluator=point, derivative=tangent) + + +def _rational_quarter_circle_cut_cell() -> CurveTrimmedPanel2D: + """Return a curved CAD-style cut cell with one rational trim edge.""" + + curved = _rational_quarter_circle_edge() + p0 = curved.point(0.0) + p1 = curved.point(1.0) + cell_corner = (0.5, 0.5) + loop = CurveLoop2D( + ( + curved, + CurveEdge2D.line(p1, cell_corner), + CurveEdge2D.line(cell_corner, p0), + ) + ) + return CurveTrimmedPanel2D(outer=loop) + + +def _curve_panel_local_residual( + panel: CurveTrimmedPanel2D, + target: Point2D, + *, + sigma: float, + order: int, +) -> float: + if order < 1: + raise ValueError("order must be positive") + if sigma <= 0.0: + raise ValueError("sigma must be positive") + + folded = folded_curve_quadrature_rule(panel, order=order) + total = 0.0 + for source, weight in zip(folded.rule.points, folded.rule.weights, strict=True): + total += ( + ewald_log_local_kernel(target, source, sigma=sigma) + * _physical_density(source) + * weight + ) + return total + + +def _quarter_circle_ray_exit_distance(target: Point2D, direction: Point2D) -> float: + """Return ray exit length for the exact rational quarter-circle cut cell.""" + + center = (0.5, 0.5) + rel = _sub(target, center) + circle_dot = _dot(rel, direction) + if circle_dot >= 0.0: + return 0.0 + + exit_distance = -2.0 * circle_dot + dx, dy = direction + if dx > 0.0: + exit_distance = min(exit_distance, (0.5 - target[0]) / dx) + if dy > 0.0: + exit_distance = min(exit_distance, (0.5 - target[1]) / dy) + return max(0.0, exit_distance) + + +def _quarter_circle_exact_boundary_local_residual( + target: Point2D, + *, + sigma: float, + angular_order: int, +) -> float: + """Integrate the local residual using exact target-centered ray clipping.""" + + if sigma <= 0.0: + raise ValueError("sigma must be positive") + if angular_order < 1: + raise ValueError("angular_order must be positive") + + p0 = (0.0, 0.5) + p1 = (0.5, 0.0) + center = (0.5, 0.5) + tangent_low = -0.25 * pi + tangent_high = 0.75 * pi + intervals = ( + ( + tangent_low, + atan2(p1[1] - target[1], p1[0] - target[0]), + "circle_start", + ), + ( + atan2(p1[1] - target[1], p1[0] - target[0]), + atan2(center[1] - target[1], center[0] - target[0]), + "right", + ), + ( + atan2(center[1] - target[1], center[0] - target[0]), + atan2(p0[1] - target[1], p0[0] - target[0]), + "top", + ), + ( + atan2(p0[1] - target[1], p0[0] - target[0]), + tangent_high, + "circle_end", + ), + ) + total = 0.0 + for theta_start, theta_end, exit_kind in intervals: + + def integrand(theta: float, *, kind: str = exit_kind) -> float: + return _quarter_circle_ray_integrand( + target, + theta, + sigma=sigma, + exit_kind=kind, + ) + + total += _adaptive_simpson( + integrand, + theta_start, + theta_end, + abs_tol=1.0e-16, + max_depth=max(20, angular_order), + ) + return total + + +def _quarter_circle_ray_integrand( + target: Point2D, + theta: float, + *, + sigma: float, + exit_kind: str, +) -> float: + direction = (cos(theta), sin(theta)) + if exit_kind in ("circle_start", "circle_end"): + exit_distance = _quarter_circle_ray_exit_distance(target, direction) + elif exit_kind == "right": + exit_distance = (0.5 - target[0]) / direction[0] + else: + exit_distance = (0.5 - target[1]) / direction[1] + if exit_distance <= 0.0: + return 0.0 + density_constant = _physical_density(target) + density_linear = ( + 0.25 * direction[0] + - 0.15 * direction[1] + + 0.1 * (target[0] * direction[1] + target[1] * direction[0]) + ) + density_quadratic = 0.1 * direction[0] * direction[1] + return ( + density_constant + * _ewald_log_local_finite_radial_moment_2d( + sigma=sigma, + monomial_degree=0, + radius=exit_distance, + ) + + density_linear + * _ewald_log_local_finite_radial_moment_2d( + sigma=sigma, + monomial_degree=1, + radius=exit_distance, + ) + + density_quadratic + * _ewald_log_local_finite_radial_moment_2d( + sigma=sigma, + monomial_degree=2, + radius=exit_distance, + ) + ) + + +def _adaptive_simpson( + function: Callable[[float], float], + start: float, + end: float, + *, + abs_tol: float, + max_depth: int, +) -> float: + midpoint = 0.5 * (start + end) + f_start = function(start) + f_mid = function(midpoint) + f_end = function(end) + whole = _simpson_estimate(start, end, f_start, f_mid, f_end) + return _adaptive_simpson_recursive( + function, + start, + midpoint, + end, + f_start, + f_mid, + f_end, + whole, + abs_tol, + max_depth, + ) + + +def _simpson_estimate( + start: float, + end: float, + f_start: float, + f_mid: float, + f_end: float, +) -> float: + return (end - start) * (f_start + 4.0 * f_mid + f_end) / 6.0 + + +def _adaptive_simpson_recursive( + function: Callable[[float], float], + start: float, + midpoint: float, + end: float, + f_start: float, + f_mid: float, + f_end: float, + whole: float, + abs_tol: float, + depth_remaining: int, +) -> float: + left_mid = 0.5 * (start + midpoint) + right_mid = 0.5 * (midpoint + end) + f_left_mid = function(left_mid) + f_right_mid = function(right_mid) + left = _simpson_estimate(start, midpoint, f_start, f_left_mid, f_mid) + right = _simpson_estimate(midpoint, end, f_mid, f_right_mid, f_end) + delta = left + right - whole + if depth_remaining <= 0 or abs(delta) <= 15.0 * abs_tol: + return left + right + delta / 15.0 + return _adaptive_simpson_recursive( + function, + start, + left_mid, + midpoint, + f_start, + f_left_mid, + f_mid, + left, + 0.5 * abs_tol, + depth_remaining - 1, + ) + _adaptive_simpson_recursive( + function, + midpoint, + right_mid, + end, + f_mid, + f_right_mid, + f_end, + right, + 0.5 * abs_tol, + depth_remaining - 1, + ) + + +def _ewald_log_local_finite_radial_moment_2d( + *, + sigma: float, + monomial_degree: int, + radius: float, +) -> float: + """Return ``int_0^radius K_sigma(r) r^(monomial_degree + 1) dr``.""" + + if sigma <= 0.0: + raise ValueError("sigma must be positive") + if radius <= 0.0: + return 0.0 + if monomial_degree not in (0, 1, 2): + raise ValueError("finite radial moment supports degrees 0, 1, and 2") + + upper = (radius / sigma) ** 2 + sqrt_upper = sqrt(upper) + exp_upper = exp(-upper) if upper < 745.0 else 0.0 + e1_upper = exponential_integral_e1(upper) + if monomial_degree == 0: + lower_gamma = 1.0 - exp_upper + q_plus_one = 1.0 + elif monomial_degree == 1: + lower_gamma = 0.5 * sqrt(pi) * erf(sqrt_upper) - sqrt_upper * exp_upper + q_plus_one = 1.5 + else: + lower_gamma = 1.0 - exp_upper * (1.0 + upper) + q_plus_one = 2.0 + + integral_u = (upper**q_plus_one * e1_upper + lower_gamma) / q_plus_one + return sigma ** (monomial_degree + 2) * integral_u / (8.0 * pi) + + +def _quarter_circle_boundary_frame() -> tuple[Point2D, Point2D, Point2D, float]: + panel = _rational_quarter_circle_cut_cell() + target = panel.outer.edges[0].point(0.5) + inv_sqrt2 = 1.0 / sqrt(2.0) + tangent = (inv_sqrt2, -inv_sqrt2) + inward_normal = (inv_sqrt2, inv_sqrt2) + return target, tangent, inward_normal, 0.5 + + +def _quarter_circle_taylor_graph_coefficients( + *, radius: float, boundary_order: int +) -> tuple[float, ...]: + if boundary_order < 2 or boundary_order % 2 != 0: + raise ValueError("boundary_order must be a positive even order at least 2") + + coefficients: list[float] = [] + binomial = 1.0 + for power in range(1, boundary_order // 2 + 1): + binomial *= (1.5 - float(power)) / float(power) + coefficient = -binomial * ((-1.0) ** power) / (radius ** (2 * power - 1)) + coefficients.append(coefficient) + return tuple(coefficients) + + +def _quarter_circle_taylor_graph_value( + s_coord: float, + *, + coefficients: tuple[float, ...], +) -> float: + value = 0.0 + s_squared = s_coord * s_coord + power = s_squared + for coefficient in coefficients: + value += coefficient * power + power *= s_squared + return value + + +def _taylor_graph_exit_distance( + theta: float, + *, + coefficients: tuple[float, ...], + sigma: float, +) -> float | None: + cos_theta = cos(theta) + sin_theta = sin(theta) + if sin_theta <= 0.0: + return 0.0 + if abs(cos_theta) <= 1.0e-14: + return None + + def boundary_residual(radius: float) -> float: + return radius * sin_theta - _quarter_circle_taylor_graph_value( + radius * cos_theta, + coefficients=coefficients, + ) + + low = 0.0 + high = max(sigma, 1.0e-3) + while boundary_residual(high) > 0.0: + high *= 2.0 + if high > 8.0: + return None + + for _ in range(96): + mid = 0.5 * (low + high) + if boundary_residual(mid) > 0.0: + low = mid + else: + high = mid + return high + + +def _local_boundary_ray_integral( + target: Point2D, + direction: Point2D, + *, + sigma: float, + exit_distance: float | None, +) -> float: + density_constant = _physical_density(target) + density_linear = ( + 0.25 * direction[0] + - 0.15 * direction[1] + + 0.1 * (target[0] * direction[1] + target[1] * direction[0]) + ) + density_quadratic = 0.1 * direction[0] * direction[1] + + def radial_moment(degree: int) -> float: + if exit_distance is None: + return ewald_log_local_radial_moment_2d( + sigma=sigma, + monomial_degree=degree, + ) + return _ewald_log_local_finite_radial_moment_2d( + sigma=sigma, + monomial_degree=degree, + radius=exit_distance, + ) + + return ( + density_constant * radial_moment(0) + + density_linear * radial_moment(1) + + density_quadratic * radial_moment(2) + ) + + +def _quarter_circle_exact_local_boundary_model( + *, sigma: float, angular_order: int +) -> float: + target, tangent, inward_normal, radius = _quarter_circle_boundary_frame() + if angular_order < 1: + raise ValueError("angular_order must be positive") + + def integrand(theta: float) -> float: + direction = ( + cos(theta) * tangent[0] + sin(theta) * inward_normal[0], + cos(theta) * tangent[1] + sin(theta) * inward_normal[1], + ) + exit_distance = 2.0 * radius * sin(theta) + return _local_boundary_ray_integral( + target, + direction, + sigma=sigma, + exit_distance=exit_distance, + ) + + return _adaptive_simpson( + integrand, + 0.0, + pi, + abs_tol=1.0e-16, + max_depth=max(20, angular_order), + ) + + +def _quarter_circle_taylor_boundary_model( + *, sigma: float, boundary_order: int, angular_order: int +) -> float: + target, tangent, inward_normal, radius = _quarter_circle_boundary_frame() + coefficients = _quarter_circle_taylor_graph_coefficients( + radius=radius, + boundary_order=boundary_order, + ) + + def integrand(theta: float) -> float: + direction = ( + cos(theta) * tangent[0] + sin(theta) * inward_normal[0], + cos(theta) * tangent[1] + sin(theta) * inward_normal[1], + ) + exit_distance = _taylor_graph_exit_distance( + theta, + coefficients=coefficients, + sigma=sigma, + ) + return _local_boundary_ray_integral( + target, + direction, + sigma=sigma, + exit_distance=exit_distance, + ) + + return _adaptive_simpson( + integrand, + 0.0, + pi, + abs_tol=1.0e-16, + max_depth=max(20, angular_order), + ) + + +def run_taylor_boundary_model_experiment( + *, + sigmas: tuple[float, ...] = (0.01, 0.02), + boundary_orders: tuple[int, ...] = (2, 4, 6, 8, 10, 12), + angular_order: int = 32, +) -> TaylorBoundaryModelReport: + """Compare high-order local boundary Taylor models to exact circular geometry.""" + + samples: list[TaylorBoundaryModelSample] = [] + for sigma in sigmas: + exact_value = _quarter_circle_exact_local_boundary_model( + sigma=sigma, + angular_order=angular_order, + ) + denominator = max(abs(exact_value), 1.0e-300) + for boundary_order in boundary_orders: + taylor_value = _quarter_circle_taylor_boundary_model( + sigma=sigma, + boundary_order=boundary_order, + angular_order=angular_order, + ) + abs_error = abs(taylor_value - exact_value) + samples.append( + TaylorBoundaryModelSample( + sigma=sigma, + boundary_order=boundary_order, + exact_value=exact_value, + taylor_value=taylor_value, + abs_error=abs_error, + rel_error=abs_error / denominator, + ) + ) + return TaylorBoundaryModelReport( + sigmas=sigmas, + boundary_orders=boundary_orders, + angular_order=angular_order, + samples=tuple(samples), + ) + + +def _local_half_plane_boundary_model(*, sigma: float, angular_order: int) -> float: + target, tangent, inward_normal, _radius = _quarter_circle_boundary_frame() + + def integrand(theta: float) -> float: + direction = ( + cos(theta) * tangent[0] + sin(theta) * inward_normal[0], + cos(theta) * tangent[1] + sin(theta) * inward_normal[1], + ) + return _local_boundary_ray_integral( + target, + direction, + sigma=sigma, + exit_distance=None, + ) + + return _adaptive_simpson( + integrand, + 0.0, + pi, + abs_tol=1.0e-16, + max_depth=max(20, angular_order), + ) + + +def _graph_strip_correction( + *, + sigma: float, + boundary_order: int, + strip_order: int, +) -> float: + target, tangent, inward_normal, radius = _quarter_circle_boundary_frame() + coefficients = _quarter_circle_taylor_graph_coefficients( + radius=radius, + boundary_order=boundary_order, + ) + nodes, weights = gauss_legendre_01(strip_order) + total = 0.0 + for s_node, s_weight in zip(nodes, weights, strict=True): + s_coord = radius * (2.0 * s_node - 1.0) + ds_weight = 2.0 * radius * s_weight + graph_height = _quarter_circle_taylor_graph_value( + s_coord, + coefficients=coefficients, + ) + if graph_height <= 0.0: + continue + for n_node, n_weight in zip(nodes, weights, strict=True): + n_coord = graph_height * n_node + source = ( + target[0] + s_coord * tangent[0] + n_coord * inward_normal[0], + target[1] + s_coord * tangent[1] + n_coord * inward_normal[1], + ) + if hypot(source[0] - target[0], source[1] - target[1]) <= _COINCIDENT_TOL: + continue + total += ( + ewald_log_local_kernel(target, source, sigma=sigma) + * _physical_density(source) + * ds_weight + * graph_height + * n_weight + ) + return total + + +def _quarter_circle_graph_boundary_model( + *, + sigma: float, + boundary_order: int, + strip_order: int, + angular_order: int, +) -> float: + return _local_half_plane_boundary_model( + sigma=sigma, + angular_order=angular_order, + ) - _graph_strip_correction( + sigma=sigma, + boundary_order=boundary_order, + strip_order=strip_order, + ) + + +def run_graph_boundary_model_experiment( + *, + sigmas: tuple[float, ...] = (0.005, 0.01, 0.02), + boundary_orders: tuple[int, ...] = (4, 6, 8), + strip_orders: tuple[int, ...] = (32, 64, 96), + angular_order: int = 32, +) -> GraphBoundaryModelReport: + """Compare graph-strip boundary corrections against exact circular geometry.""" + + samples: list[GraphBoundaryModelSample] = [] + for sigma in sigmas: + exact_value = _quarter_circle_exact_local_boundary_model( + sigma=sigma, + angular_order=angular_order, + ) + denominator = max(abs(exact_value), 1.0e-300) + for boundary_order in boundary_orders: + for strip_order in strip_orders: + graph_value = _quarter_circle_graph_boundary_model( + sigma=sigma, + boundary_order=boundary_order, + strip_order=strip_order, + angular_order=angular_order, + ) + abs_error = abs(graph_value - exact_value) + samples.append( + GraphBoundaryModelSample( + sigma=sigma, + boundary_order=boundary_order, + strip_order=strip_order, + exact_value=exact_value, + graph_value=graph_value, + abs_error=abs_error, + rel_error=abs_error / denominator, + ) + ) + return GraphBoundaryModelReport( + sigmas=sigmas, + boundary_orders=boundary_orders, + strip_orders=strip_orders, + angular_order=angular_order, + samples=tuple(samples), + ) + + +def _lagrange_interpolate( + nodes: tuple[float, ...], values: tuple[float, ...], point: float +) -> float: + if len(nodes) != len(values): + raise ValueError("nodes and values must have the same length") + total = 0.0 + for i, node_i in enumerate(nodes): + basis = 1.0 + for j, node_j in enumerate(nodes): + if i == j: + continue + basis *= (point - node_j) / (node_i - node_j) + total += values[i] * basis + return total + + +def run_precomputed_boundary_table_experiment( + *, + table_sigmas: tuple[float, ...] = ( + 0.004, + 0.008, + 0.012, + 0.016, + 0.020, + 0.024, + 0.028, + ), + eval_sigmas: tuple[float, ...] = (0.010, 0.018, 0.026), + boundary_order: int = 8, + angular_order: int = 32, +) -> PrecomputedBoundaryTableReport: + """Interpolate precomputed high-order smooth-boundary moment values.""" + + table_values = tuple( + _quarter_circle_taylor_boundary_model( + sigma=sigma, + boundary_order=boundary_order, + angular_order=angular_order, + ) + for sigma in table_sigmas + ) + samples: list[PrecomputedBoundaryTableSample] = [] + for sigma in eval_sigmas: + exact_value = _quarter_circle_taylor_boundary_model( + sigma=sigma, + boundary_order=boundary_order, + angular_order=angular_order, + ) + table_value = _lagrange_interpolate(table_sigmas, table_values, sigma) + abs_error = abs(table_value - exact_value) + samples.append( + PrecomputedBoundaryTableSample( + sigma=sigma, + exact_value=exact_value, + table_value=table_value, + abs_error=abs_error, + rel_error=abs_error / max(abs(exact_value), 1.0e-300), + ) + ) + return PrecomputedBoundaryTableReport( + table_sigmas=table_sigmas, + eval_sigmas=eval_sigmas, + boundary_order=boundary_order, + angular_order=angular_order, + samples=tuple(samples), + ) + + +def run_curved_boundary_local_model_experiment( + *, + sigmas: tuple[float, ...] = (0.02, 0.04, 0.08), + orders: tuple[int, ...] = (12, 20, 32), + reference_order: int = 96, +) -> CurvedBoundaryLocalModelReport: + """Test local residual models on an exact curved rational cut cell.""" + + panel = _rational_quarter_circle_cut_cell() + trim = panel.outer.edges[0] + target = trim.point(0.5) + + samples: list[CurvedBoundaryLocalModelSample] = [] + for sigma in sigmas: + reference = _curve_panel_local_residual( + panel, + target, + sigma=sigma, + order=reference_order, + ) + curved_boundary = _quarter_circle_exact_boundary_local_residual( + target, + sigma=sigma, + angular_order=reference_order, + ) + density_at_target = _physical_density(target) + half_plane = ( + ewald_log_local_sector_monomial_moment_2d( + sigma=sigma, + x_power=0, + y_power=0, + theta_start=0.0, + theta_end=pi, + ) + * density_at_target + ) + full_plane = ( + ewald_log_local_sector_monomial_moment_2d( + sigma=sigma, + x_power=0, + y_power=0, + theta_start=0.0, + theta_end=2.0 * pi, + ) + * density_at_target + ) + denominator = max(abs(reference), 1.0e-300) + for order in orders: + folded_value = _curve_panel_local_residual( + panel, + target, + sigma=sigma, + order=order, + ) + folded_error = abs(folded_value - reference) + curved_error = abs(curved_boundary - reference) + half_error = abs(half_plane - reference) + full_error = abs(full_plane - reference) + samples.append( + CurvedBoundaryLocalModelSample( + sigma=sigma, + order=order, + reference_order=reference_order, + target=target, + folded_value=folded_value, + reference_value=reference, + curved_boundary_value=curved_boundary, + half_plane_value=half_plane, + full_plane_value=full_plane, + folded_abs_error=folded_error, + curved_boundary_abs_error=curved_error, + half_plane_abs_error=half_error, + full_plane_abs_error=full_error, + curved_boundary_rel_error=curved_error / denominator, + half_plane_rel_error=half_error / denominator, + full_plane_rel_error=full_error / denominator, + ) + ) + + return CurvedBoundaryLocalModelReport( + sigmas=sigmas, + orders=orders, + reference_order=reference_order, + samples=tuple(samples), + ) + + +def _analytic_full_space_local_residual(target: Point2D, *, sigma: float) -> float: + """Return the leading full-space Taylor moment for the 2D local residual. + + For ``G_local = E1(r^2/sigma^2)/(4*pi)``, the zeroth full-space moment is + ``sigma^2/4``. The current experiment density has zero Laplacian, so the + second-order isotropic Taylor correction vanishes. + """ + + if sigma <= 0.0: + raise ValueError("sigma must be positive") + return (sigma * sigma / 4.0) * _physical_density(target) + + +def _panel_potential( + fans: tuple[FanTemplateMap2D, ...], + target: Point2D, + *, + order: int, + sigma: float | None, + smooth_only: bool, + local_only: bool = False, +) -> float: + if smooth_only and local_only: + raise ValueError("smooth_only and local_only are mutually exclusive") + + total = 0.0 + nodes, weights = gauss_legendre_01(order) + for fan in fans: + for r, wr in zip(nodes, weights, strict=True): + jacobian = fan.signed_jacobian(r) + for t, wt in zip(nodes, weights, strict=True): + source = fan.point(r, t) + if smooth_only: + if sigma is None: + raise ValueError("sigma is required for smooth split") + kernel = ewald_log_smooth_kernel(target, source, sigma=sigma) + elif local_only: + if sigma is None: + raise ValueError("sigma is required for local split") + kernel = ewald_log_local_kernel(target, source, sigma=sigma) + else: + if ( + hypot(target[0] - source[0], target[1] - source[1]) + <= _COINCIDENT_TOL + ): + raise ValueError( + "target coincides with a source quadrature node" + ) + kernel = laplace_log_kernel(target, source) + total += kernel * _physical_density(source) * jacobian * wr * wt + return total + + +def run_dmk_split_nearfield_experiment( + *, + orders: tuple[int, ...] = (4, 6, 8, 10, 14, 18), + sigmas: tuple[float, ...] = (0.08, 0.16, 0.32), + reference_order: int = 96, +) -> DmkSplitExperimentReport: + """Run the DMK-style smooth/local split experiment on straight fan panels.""" + + fixtures: tuple[ + tuple[str, tuple[Point2D, ...], tuple[tuple[str, Point2D], ...]], ... + ] = ( + ( + "convex_quad", + ((0.0, 0.0), (1.0, 0.0), (0.95, 0.42), (0.18, 0.86)), + ( + ("interior", (0.55, 0.32)), + ("boundary_near_inside", (0.54, 0.025)), + ("boundary_near_outside", (0.54, -0.025)), + ("vertex_near_inside", (0.05, 0.045)), + ), + ), + ( + "skew_quad", + ((0.0, 0.0), (1.15, 0.02), (0.72, 0.78), (0.08, 0.62)), + ( + ("interior", (0.50, 0.34)), + ("boundary_near_inside", (0.56, 0.045)), + ("boundary_near_outside", (0.56, -0.025)), + ("vertex_near_inside", (0.06, 0.05)), + ), + ), + ) + + seed_quality: list[SeedQualitySample] = [] + samples: list[DmkSplitSample] = [] + for fixture_name, nodes, targets in fixtures: + panel = TrimmedPanel2D(outer=PanelLoop2D(nodes)) + seed_options = ( + ("interior_anchor", select_interior_anchor(panel)), + ("node_barycenter", _node_barycenter(nodes)), + ) + for seed_mode, seed in seed_options: + seed_quality.append( + _panel_seed_quality( + fixture=fixture_name, + seed_mode=seed_mode, + nodes=nodes, + seed=seed, + ) + ) + fans = _panel_fans(nodes, seed=seed) + full_references = { + target_label: _panel_potential( + fans, + target, + order=reference_order, + sigma=None, + smooth_only=False, + ) + for target_label, target in targets + } + for sigma in sigmas: + smooth_references = { + target_label: _panel_potential( + fans, + target, + order=reference_order, + sigma=sigma, + smooth_only=True, + ) + for target_label, target in targets + } + local_references = { + target_label: _panel_potential( + fans, + target, + order=reference_order, + sigma=sigma, + smooth_only=False, + local_only=True, + ) + for target_label, target in targets + } + for target_label, target in targets: + full_reference = full_references[target_label] + smooth_reference = smooth_references[target_label] + local_reference = local_references[target_label] + ratio_denominator = max(abs(full_reference), 1.0e-300) + for order in orders: + full_value = _panel_potential( + fans, + target, + order=order, + sigma=None, + smooth_only=False, + ) + smooth_value = _panel_potential( + fans, + target, + order=order, + sigma=sigma, + smooth_only=True, + ) + local_value = _panel_potential( + fans, + target, + order=order, + sigma=sigma, + smooth_only=False, + local_only=True, + ) + analytic_local_value = _analytic_full_space_local_residual( + target, + sigma=sigma, + ) + reconstructed_value = smooth_value + local_value + analytic_reconstructed_value = ( + smooth_value + analytic_local_value + ) + samples.append( + DmkSplitSample( + fixture=fixture_name, + seed_mode=seed_mode, + seed=seed, + target_label=target_label, + target=target, + sigma=sigma, + order=order, + full_value=full_value, + smooth_value=smooth_value, + local_value=local_value, + analytic_local_value=analytic_local_value, + reconstructed_value=reconstructed_value, + analytic_reconstructed_value=analytic_reconstructed_value, + full_reference=full_reference, + smooth_reference=smooth_reference, + local_reference=local_reference, + full_abs_error=abs(full_value - full_reference), + smooth_abs_error=abs(smooth_value - smooth_reference), + local_abs_error=abs(local_value - local_reference), + analytic_local_abs_error=abs( + analytic_local_value - local_reference + ), + reconstructed_abs_error=abs( + reconstructed_value - full_reference + ), + analytic_reconstructed_abs_error=abs( + analytic_reconstructed_value - full_reference + ), + local_to_full_ratio=abs(local_reference) + / ratio_denominator, + ) + ) + + return DmkSplitExperimentReport( + reference_order=reference_order, + orders=orders, + sigmas=sigmas, + seed_quality=tuple(seed_quality), + samples=tuple(samples), + ) + + def expected_scaled_laplace_point_potential( base_potential: float, source_mass: float, diff --git a/tests/evals/test_nearfield_templates.py b/tests/evals/test_nearfield_templates.py index d3c60e1..e84eeea 100644 --- a/tests/evals/test_nearfield_templates.py +++ b/tests/evals/test_nearfield_templates.py @@ -1,18 +1,70 @@ from __future__ import annotations +from math import cos, exp, log, pi, sin + import pytest from cutkit.evals import ( FanTemplateMap2D, diagonal_remainder_sample, + ewald_log_local_kernel, + ewald_log_local_radial_moment_2d, + ewald_log_local_sector_monomial_moment_2d, + ewald_log_smooth_kernel, expected_scaled_laplace_point_potential, point_target_laplace_potential, + point_target_smooth_ewald_potential, + run_curved_boundary_local_model_experiment, + run_dmk_split_nearfield_experiment, + run_graph_boundary_model_experiment, run_nearfield_template_experiment, + run_precomputed_boundary_table_experiment, + run_taylor_boundary_model_experiment, template_density_mass, ) from cutkit.quadrature import gauss_legendre_01 +def _polar_sector_local_moment_reference( + *, + sigma: float, + x_power: int, + y_power: int, + theta_start: float, + theta_end: float, + radial_order: int = 220, + angular_order: int = 96, +) -> float: + log_radius_min = log(sigma) - 32.0 + log_radius_max = log(8.0 * sigma) + radial_nodes, radial_weights = gauss_legendre_01(radial_order) + angular_nodes, angular_weights = gauss_legendre_01(angular_order) + log_radius_width = log_radius_max - log_radius_min + theta_width = theta_end - theta_start + total = 0.0 + for log_radius_node, log_radius_weight_base in zip( + radial_nodes, radial_weights, strict=True + ): + log_radius = log_radius_min + log_radius_width * log_radius_node + radius = exp(log_radius) + # The polar area factor contributes one radius and dr = radius d(log r). + radial_weight = log_radius_width * log_radius_weight_base * radius * radius + for theta_node, theta_weight_base in zip( + angular_nodes, angular_weights, strict=True + ): + theta = theta_start + theta_width * theta_node + theta_weight = theta_width * theta_weight_base + source = (radius * cos(theta), radius * sin(theta)) + total += ( + ewald_log_local_kernel((0.0, 0.0), source, sigma=sigma) + * source[0] ** x_power + * source[1] ** y_power + * radial_weight + * theta_weight + ) + return total + + def test_fan_map_area_and_metric_are_consistent() -> None: fan = FanTemplateMap2D( vertex=(0.0, 0.0), @@ -142,6 +194,279 @@ def test_point_target_path_rejects_coincident_source_node() -> None: point_target_laplace_potential(fan, target, order=3) +def test_ewald_log_split_reconstructs_log_kernel() -> None: + target = (0.3, 0.4) + source = (0.9, -0.2) + sigma = 0.25 + + split = ewald_log_smooth_kernel( + target, source, sigma=sigma + ) + ewald_log_local_kernel( + target, + source, + sigma=sigma, + ) + + from cutkit.evals import laplace_log_kernel + + assert split == pytest.approx(laplace_log_kernel(target, source), abs=1.0e-13) + + +def test_ewald_local_full_plane_moments_match_known_values() -> None: + sigma = 0.17 + + mass = ewald_log_local_sector_monomial_moment_2d( + sigma=sigma, + x_power=0, + y_power=0, + theta_start=0.0, + theta_end=2.0 * pi, + ) + xx_moment = ewald_log_local_sector_monomial_moment_2d( + sigma=sigma, + x_power=2, + y_power=0, + theta_start=0.0, + theta_end=2.0 * pi, + ) + xy_moment = ewald_log_local_sector_monomial_moment_2d( + sigma=sigma, + x_power=1, + y_power=1, + theta_start=0.0, + theta_end=2.0 * pi, + ) + + assert mass == pytest.approx(sigma * sigma / 4.0, rel=1.0e-13) + assert xx_moment == pytest.approx(sigma**4 / 16.0, rel=1.0e-13) + assert xy_moment == pytest.approx(0.0, abs=1.0e-19) + + +def test_ewald_local_half_plane_moments_match_polar_reference() -> None: + sigma = 0.13 + theta_start = 0.0 + theta_end = pi + + for x_power, y_power in ((0, 0), (1, 0), (0, 1), (2, 0), (0, 2)): + analytic = ewald_log_local_sector_monomial_moment_2d( + sigma=sigma, + x_power=x_power, + y_power=y_power, + theta_start=theta_start, + theta_end=theta_end, + ) + reference = _polar_sector_local_moment_reference( + sigma=sigma, + x_power=x_power, + y_power=y_power, + theta_start=theta_start, + theta_end=theta_end, + ) + assert analytic == pytest.approx(reference, rel=5.0e-4, abs=1.0e-12) + + mass = ewald_log_local_sector_monomial_moment_2d( + sigma=sigma, + x_power=0, + y_power=0, + theta_start=theta_start, + theta_end=theta_end, + ) + assert mass == pytest.approx(sigma * sigma / 8.0, rel=1.0e-13) + + +def test_ewald_local_wedge_moments_match_polar_reference() -> None: + sigma = 0.11 + theta_start = -0.2 + theta_end = 0.85 * pi + + for x_power, y_power in ((0, 0), (1, 0), (0, 1), (2, 0), (1, 1)): + analytic = ewald_log_local_sector_monomial_moment_2d( + sigma=sigma, + x_power=x_power, + y_power=y_power, + theta_start=theta_start, + theta_end=theta_end, + ) + reference = _polar_sector_local_moment_reference( + sigma=sigma, + x_power=x_power, + y_power=y_power, + theta_start=theta_start, + theta_end=theta_end, + ) + assert analytic == pytest.approx(reference, rel=5.0e-4, abs=1.0e-12) + + mass = ewald_log_local_sector_monomial_moment_2d( + sigma=sigma, + x_power=0, + y_power=0, + theta_start=theta_start, + theta_end=theta_end, + ) + opening_angle = theta_end - theta_start + assert mass == pytest.approx(opening_angle * sigma * sigma / (8.0 * pi)) + + +def test_ewald_local_radial_moment_rejects_invalid_inputs() -> None: + with pytest.raises(ValueError, match="sigma must be positive"): + ewald_log_local_radial_moment_2d(sigma=0.0, monomial_degree=0) + with pytest.raises(ValueError, match="monomial_degree must be nonnegative"): + ewald_log_local_radial_moment_2d(sigma=0.1, monomial_degree=-1) + + +def test_curved_boundary_local_model_beats_full_plane_on_cad_style_cell() -> None: + report = run_curved_boundary_local_model_experiment( + sigmas=(0.02, 0.04), + orders=(16, 32), + reference_order=80, + ) + + finest_by_sigma = { + sample.sigma: sample for sample in report.samples if sample.order == 32 + } + + assert len(finest_by_sigma) == 2 + for sample in finest_by_sigma.values(): + assert sample.reference_value > 0.0 + assert sample.half_plane_abs_error < sample.full_plane_abs_error + assert sample.half_plane_rel_error < 0.03 + assert sample.full_plane_rel_error > 1.0 + + +def test_curved_boundary_exact_model_reaches_machine_precision() -> None: + coarse = run_curved_boundary_local_model_experiment( + sigmas=(0.02, 0.04), + orders=(16,), + reference_order=24, + ) + fine = run_curved_boundary_local_model_experiment( + sigmas=(0.02, 0.04), + orders=(16,), + reference_order=32, + ) + + coarse_by_sigma = {sample.sigma: sample for sample in coarse.samples} + fine_by_sigma = {sample.sigma: sample for sample in fine.samples} + + for sigma, fine_sample in fine_by_sigma.items(): + coarse_value = coarse_by_sigma[sigma].curved_boundary_value + assert fine_sample.curved_boundary_value == pytest.approx( + coarse_value, + rel=1.0e-10, + abs=1.0e-14, + ) + + +def test_high_order_taylor_boundary_model_reaches_machine_precision() -> None: + report = run_taylor_boundary_model_experiment( + sigmas=(0.02,), + boundary_orders=(2, 4, 6, 8), + angular_order=32, + ) + errors = {sample.boundary_order: sample.rel_error for sample in report.samples} + + assert errors[4] < errors[2] + assert errors[6] < errors[4] + assert errors[8] < 1.0e-10 + + +def test_high_order_taylor_boundary_model_handles_smaller_windows() -> None: + report = run_taylor_boundary_model_experiment( + sigmas=(0.005, 0.01), + boundary_orders=(4, 6), + angular_order=32, + ) + + finest_errors = { + sample.sigma: sample.rel_error + for sample in report.samples + if sample.boundary_order == 6 + } + + assert finest_errors[0.005] < 1.0e-10 + assert finest_errors[0.01] < 1.0e-10 + + +def test_precomputed_boundary_table_interpolates_to_machine_precision() -> None: + report = run_precomputed_boundary_table_experiment( + table_sigmas=(0.004, 0.008, 0.012, 0.016, 0.020, 0.024, 0.028), + eval_sigmas=(0.010, 0.018, 0.026), + boundary_order=8, + angular_order=32, + ) + + assert max(sample.rel_error for sample in report.samples) < 1.0e-10 + + +def test_direct_graph_strip_quadrature_is_not_enough_for_machine_precision() -> None: + report = run_graph_boundary_model_experiment( + sigmas=(0.02,), + boundary_orders=(8,), + strip_orders=(64, 128), + angular_order=32, + ) + + errors = {sample.strip_order: sample.rel_error for sample in report.samples} + + assert errors[128] < errors[64] + assert errors[128] > 1.0e-10 + + +def test_curved_boundary_folded_reference_converges_with_order() -> None: + report = run_curved_boundary_local_model_experiment( + sigmas=(0.04,), + orders=(10, 16, 32), + reference_order=80, + ) + + errors = {sample.order: sample.folded_abs_error for sample in report.samples} + + assert errors[32] < errors[16] + assert errors[16] < errors[10] + + +def test_smoothed_ewald_potential_allows_coincident_quadrature_node() -> None: + fan = FanTemplateMap2D( + vertex=(0.0, 0.0), + edge_start=(1.0, 0.0), + edge_end=(0.35, 0.9), + ) + nodes, _weights = gauss_legendre_01(3) + target = fan.point(nodes[1], nodes[1]) + + value = point_target_smooth_ewald_potential( + fan, + target, + order=3, + sigma=0.2, + ) + + assert value == pytest.approx(value) + + +def test_dmk_split_experiment_records_local_residual_reconstruction() -> None: + report = run_dmk_split_nearfield_experiment( + orders=(4,), + sigmas=(0.16,), + reference_order=8, + ) + + sample = report.samples[0] + assert sample.reconstructed_value == pytest.approx( + sample.smooth_value + sample.local_value + ) + assert sample.analytic_reconstructed_value == pytest.approx( + sample.smooth_value + sample.analytic_local_value + ) + assert sample.local_reference == pytest.approx( + sample.full_reference - sample.smooth_reference, + abs=1.0e-12, + ) + assert sample.reconstructed_abs_error >= 0.0 + assert sample.analytic_local_abs_error >= 0.0 + assert sample.analytic_reconstructed_abs_error >= 0.0 + + def test_baseline_experiment_records_feasibility_checks() -> None: result = run_nearfield_template_experiment(order=5, scale_factor=1.5) From fb11af8249007f16cda045b4fafc91e173e99efc Mon Sep 17 00:00:00 2001 From: xywei Date: Mon, 15 Jun 2026 19:33:42 -0500 Subject: [PATCH 15/15] docs: add complement QBX near-boundary route Record the QBX-style fix for complement-subtraction singular quadrature: when the target is too close to the complement boundary for efficient direct quadrature, compute a smooth complement expansion from a separated center and evaluate it at the near-boundary target.\n\nThis keeps the self singularity in the full-box table while avoiding nearly singular folded quadrature over the complement.\n\nTested with uv sync --extra dev uv run prek install --overwrite --hook-type pre-commit --hook-type pre-push Overwriting existing hook at `/home/xywei/Work/fmm/cutkit/.git/hooks/pre-commit` prek installed at `/home/xywei/Work/fmm/cutkit/.git/hooks/pre-commit` Overwriting existing hook at `/home/xywei/Work/fmm/cutkit/.git/hooks/pre-push` prek installed at `/home/xywei/Work/fmm/cutkit/.git/hooks/pre-push` uv run ruff format --check . 88 files already formatted uv run ruff check . All checks passed! uv run mypy src tests Success: no issues found in 75 source files uv run python scripts/check_architecture.py Architecture check passed. uv run python scripts/check_docs_freshness.py Docs freshness check passed. uv run pytest -q ........................................................................ [ 18%] ........................................................................ [ 36%] ........................................................................ [ 54%] .................................................................s...... [ 72%] ........................................................................ [ 90%] ...................................... [100%] uv run python scripts/run_cutpanel_eval.py name | area | cut_fraction | area_err | max_moment_err | status --- | --- | --- | --- | --- | --- unit-square | 1.000000 | 1.000000 | 6.661e-16 | 1.110e-15 | PASS square-with-hole | 0.750000 | 0.750000 | 5.551e-16 | 3.886e-16 | PASS rect-with-thin-hole | 1.980000 | 0.990000 | 1.554e-15 | 3.109e-15 | PASS square-with-seam-adjacent-slot | 0.998200 | 0.998200 | 6.661e-16 | 1.998e-15 | PASS square-with-ultra-thin-frame | 0.011964 | 0.011964 | 1.076e-16 | 5.065e-16 | PASS imported-prod-courtyard-double-cutout | 12.040000 | 0.802667 | 1.066e-14 | 9.237e-14 | PASS imported-prod-long-corridor-offset-room | 11.620000 | 0.806944 | 8.882e-15 | 3.411e-13 | PASS imported-prod-staggered-triple-core | 23.980000 | 0.749375 | 1.776e-14 | 1.023e-12 | PASS imported-prod-offset-courtyard-barbell | 18.360000 | 0.690226 | 1.776e-14 | 3.411e-13 | PASS imported-prod-service-spine-notches | 21.780000 | 0.806667 | 1.066e-14 | 1.933e-12 | PASS fuzz-candidate-dual-slim-slots-a | 0.894400 | 0.894400 | 5.551e-16 | 1.443e-15 | PASS fuzz-candidate-near-seam-pair | 0.958900 | 0.958900 | 5.551e-16 | 6.661e-16 | PASS fuzz-candidate-quad-pocket-micro | 2.686200 | 0.839438 | 1.332e-15 | 5.773e-15 | PASS fuzz-candidate-staggered-slot-grid | 3.445000 | 0.768973 | 3.553e-15 | 7.994e-15 | PASS fuzz-candidate-thin-z-corridor | 2.208000 | 0.766667 | 8.882e-16 | 3.997e-15 | PASS fuzz-candidate-triple-pocket-offset | 2.607000 | 0.846429 | 2.220e-15 | 3.553e-15 | PASS Cut-panel evaluation passed.. Entire-Checkpoint: 1db770c7d7f6 --- docs/nearfield-local-model-catalogue.md | 41 +++++++++++++++++++++---- docs/nearfield-template-experiments.md | 36 +++++++++++++++++++++- 2 files changed, 70 insertions(+), 7 deletions(-) diff --git a/docs/nearfield-local-model-catalogue.md b/docs/nearfield-local-model-catalogue.md index f667b78..c70504d 100644 --- a/docs/nearfield-local-model-catalogue.md +++ b/docs/nearfield-local-model-catalogue.md @@ -506,15 +506,43 @@ target lies on the physical boundary, when the complement closure contains the target, or when the target-complement distance is small relative to the requested accuracy and quadrature order. +For the small-distance case, the preferred fix is a QBX-style expansion of the +complement potential rather than direct high-resolution quadrature at the target. +Although the integrand `K(x,y)` is nearly singular when `x` approaches the cut +boundary, the complement potential is smooth on the physical side away from the +complement sources. Choose an expansion center `c` in `Omega_B` and a radius `r` +such that + +$$ +|x-c| < r < \operatorname{dist}(c, B\setminus\Omega). +$$ + +Then compute local expansion coefficients for + +$$ +u_{\mathrm{comp}}(x)=\int_{B\setminus\Omega}K(x,y)p_\alpha(y)\,dy +$$ + +by folded quadrature over the complement with target/center `c`, where the +coefficient integrands are smooth. The expansion is then evaluated at the actual +near-boundary target `x` and subtracted from the full-box singular moment. + +This is a volume-potential analogue of QBX: the expansion is not used to represent +the singular self term on `Omega_B`; that part is already in the full-box table. +It is used only to avoid nearly singular evaluation of the smooth complement +field near the boundary. + Routing for this route is: 1. If the target is target-separated from the complement, use full-box table minus folded complement. -2. If the target lies on one smooth boundary feature, use smooth-boundary moments - or a boundary table instead. -3. If the target lies on an edge, corner, or vertex, route to the corresponding +2. If the target is close to the complement boundary but covered by a safe + expansion ball, use full-box table minus complement-QBX evaluation. +3. If no safe complement expansion center exists and the target lies on one + smooth boundary feature, use smooth-boundary moments or a boundary table. +4. If the target lies on an edge, corner, or vertex, route to the corresponding wedge/cone model or refine. -4. If multiple complement pieces approach the target, shrink the box/window, +5. If multiple complement pieces approach the target, shrink the box/window, subdivide the geometry, or use a target-centered fallback. The complement folded quadrature should use open quadrature rules on fan or cone @@ -527,8 +555,9 @@ the smooth interior of each complement piece. This route is attractive for List 1 matrix generation because the singular part is geometry independent and can reuse the same full-box tables as uncut boxes. The geometry-dependent object is the complement moment matrix, which is smooth -under the separation condition and can be built by folded quadrature or compressed -over cut-geometry parameters. +under the separation condition and can be built by folded quadrature, by +QBX-style local expansion coefficients, or compressed over cut-geometry +parameters. ## Precomputed Smooth-Boundary Moment Scheme diff --git a/docs/nearfield-template-experiments.md b/docs/nearfield-template-experiments.md index adc1c3c..bdd22b3 100644 --- a/docs/nearfield-template-experiments.md +++ b/docs/nearfield-template-experiments.md @@ -421,6 +421,34 @@ quadrature to be efficient. Those cases still need boundary-aware moment tables, target-centered Duffy fallback, smaller boxes, or a residual window that sees only one feature. +There is a useful intermediate case before falling back to singular boundary +moments. The complement contribution + +$$ +u_{\mathrm{comp}}(x)=\int_{B\setminus\Omega}K(x,y)p_j(y)\,dy +$$ + +is smooth as a function of `x` throughout the physical cut cell, even when `x` is +very close to the cut boundary and direct quadrature of the nearly singular +integrand is expensive. This is the same opportunity used by QBX: choose an +expansion center `c` inside `Omega_B`, separated from the complement by a safe +radius, compute a local expansion of `u_comp` about `c`, and evaluate that +expansion at near-boundary target nodes covered by the expansion ball. + +In this version, the complement term becomes: + +1. folded quadrature of expansion coefficients at a center with safe separation + from `B \setminus Omega`; +2. local Taylor, harmonic, or kernel-specific QBX evaluation at the physical + target nodes; +3. subtraction from the full-box singular table. + +This preserves the main advantage of the complement route: the singularity is +still handled by the full-box table, while the geometry-dependent complement is +handled through smooth coefficient integrals. The expansion center and order must +be chosen so the expansion ball covers the target but does not cross the +complement source region. + Target-centered Duffy refolding is useful as a diagnostic or single-target fallback, but it should not be treated as the main reusable strategy. Folded decomposition does not require a fixed seed; for one point target `x` inside the @@ -573,7 +601,10 @@ models. Its expected sweet spot is a cut-box target that is interior to `Omega cap B` but close enough to the cut boundary that direct folded quadrature of the singular kernel is poor. The full-box singular table supplies the singular moment exactly for the ambient basis, and the complement integral should converge -as a smooth folded integral as long as the complement is target-separated. +as a smooth folded integral as long as the complement is target-separated. For +near-boundary targets where direct complement quadrature is nearly singular, test +the complement-QBX variant: compute expansion coefficients from folded complement +quadrature at a separated center and evaluate the expansion at the target. For each geometry and interaction case, record: @@ -600,6 +631,9 @@ For each geometry and interaction case, record: - whether a full-box singular moment minus a smooth folded complement integral is more accurate or cheaper than direct boundary-aware residual tables for target-separated cut-box interactions; +- whether complement-QBX expansions remove the close-to-boundary resolution + burden in the complement integral, and how expansion center placement/order + should be chosen relative to the cut boundary; - whether metric-only organization is sufficient for any subfamily, or whether higher-order jet coefficients dominate the correction size. - how often source-box and neighbor-box target nodes actually require the