Repository navigation
fix bugs in 3D point search / intersections - #384
Merged
jacobmerson merged 13 commits intoOct 11, 2026
Merged
jacobmerson merged 13 commits into
jacobmerson merged 13 commits into
Conversation
The grid cell -> element candidate map was built cell-major: for every grid cell, test every element. With count_and_fill_crs running the functor twice, that is O(2 * num_cells * nelems) -- 22.9e9 tet-vs-bbox tests just to build the two 10^3 grids for a 5.7M element mesh, which measured at 175.8 s, and 91.5e9 for the 20^3 grid the conservative projection uses. Invert the iteration: each element visits only the cells its own bounding box spans, which is 1-8 cells rather than all of them, then transpose the element -> cell map into the cell -> element map the search needs. Building the same two 10^3 grids now takes 6.7 s. The element's cell range is widened by one cell per axis before applying the overlap predicate. The predicate treats a shared boundary plane as an intersection, so a cell that only touches the bounding box is a candidate, and on a structured mesh whose vertices land on grid lines that is the common case, not an edge case. One cell is sufficient because a cell two or more away cannot touch the box. Rows are sorted after the atomic scatter to keep ascending element IDs. The query returns the first candidate containing the point, so row order decides which element wins when a point sits on a shared face, and the out-of-bounds path reports the row's first entry; unordered rows would make both vary run to run. Kokkos::transpose_crs is not used for the transpose: it sizes the output row count from in.numRows(), which over-allocates when cells are fewer than elements and writes out of bounds when they are more. The previous cell-major build is kept as construct_intersection_map_reference_* and used as a test oracle, asserting the fast build reproduces it exactly -- same row offsets and same entries in the same order -- in 2D and 3D at grids coarser and finer than the mesh. That test is what caught the missing boundary-touching candidates above.
The background grid resolution was hard coded and independent of the mesh: 10^3 in the Lagrange evaluator factory and 20^3 in the mesh intersection search. A point query linearly scans its cell's candidate list, so cell occupancy sets the query cost, and a fixed resolution makes that cost grow with the element count. 8000 cells over a 5.7M element mesh is ~715 elements per cell on average, and far worse in occupied cells because a stellarator fills only a fraction of its bounding box. Scale the divisions as nelems^(1/dim) so occupancy stays roughly constant, with a floor of 10 per axis for small meshes. This is only affordable now that the candidate map is built element-major; the previous cell-major build was O(num_cells * nelems), so raising the resolution made construction proportionally worse. Kokkos::transpose_crs would also have been out of bounds here, since the cell count now exceeds the element count. Behavior is unchanged for points inside the mesh. Out-of-bounds queries report the nearest candidate from a different cell, so the specific element ID they return can change.
The adjacency search walks the source mesh's dual graph, so most elements it examines share a face with one that already overlaps the target and contribute no volume. Each was rejected by running a full r3d clip and testing the resulting measure against a tolerance, in both the counting and the filling pass. On a 5.7M element same-mesh projection that is about 4 rejections per target element per pass, and rejections are roughly 80% of all clips. Add simplices_have_degenerate_overlap, which uses the face planes of both simplices as candidate separating axes: if every vertex of one lies in the closed outer half-space of some face plane of the other, the intersection is contained in that plane and has no volume. Face-adjacent elements are separated exactly by their shared face's plane, so they are settled in one or two plane tests rather than a clip. The test is deliberately conservative in two ways. It rejects only on provable separation, making no attempt to reproduce the volume tolerance the clip applies -- a pair with a small but admissible overlap has a vertex strictly inside, falls through to the clip, and gets the same verdict as before. And face planes alone are an incomplete set of separating axes, so a genuinely disjoint pair may go undetected, which only costs a clip that would have happened anyway. intersectTargets and adjBasedIntersectSearch take a use_prefilter flag purely so the tests can compare against the unfiltered path, which serves as the oracle: the resulting intersection map must be identical, offsets and indices and order, over 2D and 3D at identical, nested and non-nested resolutions. A second test pins down that the predicate actually rejects face-adjacent elements and does not reject overlapping ones, since a filter that never fired would satisfy the equivalence test while saving nothing.
When a source element shares a face plane with the target (always the case on the same mesh), r3d sees signed distances of order 1e-13 with mixed signs on that plane and splices new vertices at O(1) fractions along the edges between them. The clipped polytope keeps the exact volume because r3d integrates signed pieces, but its face graph is folded and some faces come back mirrored. Taking |volume| of each star-decomposition piece broke that cancellation and over-counted the element by up to 30%, which on the PUMI-Tally same-mesh P0->P1 case showed up as a 1.4e-6 conservation error concentrated on 118 elements. Factor the star decomposition into ForEachPolytopeStarTet, accumulate signed piece volumes, and fix the overall orientation from their total, mirroring r3d::integrate. The RHS integrator already uses the measure as a quadrature weight, so mirrored pieces now cancel as they should. Add tests: a pinned folded polytope captured from the tally mesh (the fold depends on floating-point contraction, so it cannot be regenerated reliably from coordinates), a self-clip volume check on the offending tets plus random tets at realistic scale, and same-mesh P0->P1 moment tests: lumped equals the element-nodal average to roundoff; consistent preserves zeroth and first moments.
Building a conservative projection located the target centroids with a freshly constructed grid search over the source mesh, once per pass of the two-pass intersection, while the source function space already owned an identical search for point evaluation. On the 5.7M-tet PUMI-Tally mesh that is two extra candidate-map builds and two centroid locations per operator. intersectTargets gains an overload that borrows a caller-owned GridPointSearchVariant, and both overloads now locate the centroids once and hand the seeds to the count and fill passes. The Omega_h Lagrange evaluator factory exposes its search, LagrangeFunctionSpace exposes its factory, and the RHS integrator's space constructor threads the source space's search through; the layout constructor keeps building its own. GetOwningElementIds becomes const, which it always was in effect. Tests: the map is the same set per target whether the centroids were located by the default grid or a caller-supplied coarse one (row order may differ when a centroid lies on a shared source face), a search of the wrong dimension is rejected, and the space and layout constructors of the RHS integrator produce identical integration points.
The 2D and 3D sub-simplex decompositions of a clipped r3d polytope were two unrelated code paths. 3D walked r3d's face graph and, since the last commit, accumulates signed pieces so folded faces cancel. 2D discarded the graph instead: it merged vertices within an absolute 1e-12, re-sorted the survivors by angle about the centroid, rebuilt the links, flipped the polygon if inverted, and only then fanned it. That happened to be immune to the fold, but only because angular order determines a convex polygon's boundary, which has no 3D analogue, and it ran atan2 per vertex per clip inside the kernel. Both dimensions now go through one ForEachPolytopeStarSimplex: centroid apex, signed piece measures via Omega_h::simplex_size_from_basis, orientation read off the signed total, one epsilon rule. The only dimension-specific step left is ForEachPolytopeBoundarySimplex, which yields the boundary (Dim-1)-simplices as r3d's graph encodes them: every edge of the cycle in 2D, the existing face walk fanned into triangles in 3D. That split mirrors r3d itself, whose reduce is written once per dimension. RemoveDuplicateVerticesAndFixLinks and the two per-dimension drivers are gone. Tests now run the same checks in both dimensions: clean-simplex pieces, a hand-built folded polygon alongside the pinned folded polytope, self-clip tiling on random and scaled-box meshes, same-mesh lumped P0->P1 against the element-nodal average, and consistent-mass moment preservation. The 3D results on the PUMI-Tally case are bit-identical to the previous commit.
The star decomposition lifted every boundary facet to the polytope's centroid, so a clean simplex became Dim+1 pieces and each piece carried its own quadrature rule. Nothing in the signed-sum argument needs the centroid: by the divergence theorem the signed pieces from any apex sum to the enclosed measure, folds included. Anchoring at the first vertex makes the facets incident to it vanish (the epsilon filter drops them) and leaves one piece for a clean simplex, so a same-mesh element gets one quadrature rule instead of four in 3D and three in 2D. On the 5.7M-tet PUMI-Tally case the lumped projection build drops from 81.8 s to 50.0 s and the apply from 1.6 s to 0.6 s, with the integration point set shrinking by about 4x; conservation and the nodal-average check stay at 1e-13. Results move at roundoff level in both dimensions, which the same-mesh tests cover; the clean-simplex test now pins the one-piece property.
The search routines seem to fail in 3D when some entities are overlapping. We have added some regression tests that are currently failing.
jacobmerson
marked this pull request as ready for review
September 27, 2026 07:22
jacobmerson
requested review from
abhiyanpaudel
and
a balanced review from Copilot
September 27, 2026 07:22
Collaborator
Author
|
This fixes a bunch of issues with the the search in 3D and adds a lumped mass smoother...probably would be better as two separate PRs and I will consider splitting it up if I get a chance. |
Contributor
There was a problem hiding this comment.
Copilot review overview
🟡 Changes recommended
Point-search ID encoding can return false owners, and rejected smoother construction leaks a PETSc vector.
Review effort: Balanced
Findings: 1
Open (2)
What changed in this PR
This PR improves point-search and mesh-intersection correctness/performance, introduces reusable intersection quadrature, and adds conservative mass smoothing.
Changes:
- Builds point-search candidate maps element-major and improves 3D shared-face localization.
- Adds intersection prefiltering, signed polytope decomposition, and reusable intersections/quadrature.
- Adds a mass smoother, Python bindings, and extensive regression tests.
| File | Description |
|---|---|
test/test_point_search.cpp |
Tests candidate maps and shared-face searches. |
test/test_omega_h_same_mesh_p0_to_p1.cpp |
Tests signed decomposition and same-mesh projection. |
test/test_omega_h_intersection_rhs_integrator.cpp |
Compares integrator constructors. |
test/test_omega_h_intersection_coverage.cpp |
Tests complete target coverage. |
test/test_omega_h_form_integrator_utils.cpp |
Updates polytope-boundary tests. |
test/test_mesh_intersection_prefilter.cpp |
Tests prefiltering and search reuse. |
test/test_mesh_intersection_field_transfer.cpp |
Tests shared intersections and quadrature. |
test/test_mass_smoother.cpp |
Covers mass-smoother behavior and validation. |
test/CMakeLists.txt |
Registers new tests. |
src/pcms/utility/uniform_grid.h |
Centralizes indexing and adaptive grid sizing. |
src/pcms/transfer/omega_h_intersection_rhs_integrator.hpp |
Exposes reusable quadrature APIs. |
src/pcms/transfer/omega_h_intersection_rhs_integrator.cpp |
Delegates storage to shared quadrature. |
src/pcms/transfer/omega_h_intersection_quadrature.hpp |
Defines reusable intersection quadrature. |
src/pcms/transfer/omega_h_intersection_quadrature.cpp |
Builds quadrature and source evaluators. |
src/pcms/transfer/omega_h_form_integrator_utils.hpp |
Implements signed star decomposition. |
src/pcms/transfer/omega_h_conservative_projection.hpp |
Adds reusable intersection inputs. |
src/pcms/transfer/omega_h_conservative_projection.cpp |
Reuses quadrature and direct evaluators. |
src/pcms/transfer/mesh_intersection.hpp |
Adds prefilter and intersection abstractions. |
src/pcms/transfer/mesh_intersection.cpp |
Reuses searches and applies prefiltering. |
src/pcms/transfer/mass_smoother.hpp |
Defines the mass smoother API. |
src/pcms/transfer/mass_smoother.cpp |
Implements conservative smoothing. |
src/pcms/transfer/CMakeLists.txt |
Builds new transfer components. |
src/pcms/transfer/bilinear_form_integrator.hpp |
Clarifies matrix indexing. |
src/pcms/pythonapi/bind_transfer_field.cpp |
Exposes new transfer APIs to Python. |
src/pcms/localization/point_search.h |
Extends point-search APIs and contracts. |
src/pcms/localization/point_search.cpp |
Implements element-major maps and 3D localization. |
src/pcms/field/function_space/lagrange.h |
Exposes the evaluator factory. |
src/pcms/field/function_space/lagrange.cpp |
Implements factory access. |
src/pcms/field/field_operator.hpp |
Adds an in-place operator interface. |
src/pcms/field/evaluation_request.h |
Supports known-element evaluation. |
src/pcms/field/CMakeLists.txt |
Installs the operator header. |
💡 Add a code-review agent skill or configure MCP servers for context-aware, tailored reviews. Learn more in the docs.
jacobmerson
changed the title
Add NO_INTERSECT and functions to address outside points. This fixes the problem with -0 == 0 and adds an enum to capture the case when no candidate cells exist.
jacobmerson
force-pushed
the
point-search-element-major-candidate-map
branch
from
October 11, 2026 18:59
040957a to
a896455
Compare
jacobmerson
changed the title
jacobmerson
changed the title
jacobmerson
commented
Oct 11, 2026
Co-authored-by: Jacob Merson <jacob@jacobmerson.com>
Co-authored-by: Jacob Merson <jacob@jacobmerson.com>
This branch was successfully deployed
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
This PR resolves some performance and correctness bugs in the 3D point search and intersections code.