Skip to content

Implement first set of fast approximation GWN methods#1791

Merged
jcs15c merged 40 commits intodevelopfrom
feature/spainhour/fast_gwn_triangles
Mar 23, 2026
Merged

Implement first set of fast approximation GWN methods#1791
jcs15c merged 40 commits intodevelopfrom
feature/spainhour/fast_gwn_triangles

Conversation

@jcs15c
Copy link
Copy Markdown
Contributor

@jcs15c jcs15c commented Feb 2, 2026

Summary

  • This PR enhances the GWN capabilities of quest_winding_number_2d_ex (previously quest_winding_number_ex) and quest_winding_number_3d_ex by adding options for fast approximations of the GWN field. These methods are currently implemented for collections of 3D triangles and 2D lines as described in this paper.
  • To support this, the quest_winding_number_3d_ex example now has command line arguments to triangulate the input STEP file and evaluate its GWN field, and quest_winding_number_2d_ex now has command line arguments to linearize input curves. Subsequent containment decisions will faithfully reproduce the boundary of the discrete shape, but this will not exactly match the boundary of the curved surface. This leads to classification errors for points close to the surface.
  • New tests for 2D and 3D fast-approximate GWN methods have been added.

Algorithm Details

  • Briefly, these fast approximations work by adding each piece of the shape to a BVH spatial index, which is traversed at the given query point. During traversal:
    • If the query point is far from the bounding box of an internal node, the contribution of all its children to the GWN is be approximated rapidly using a Taylor expansion.
    • If the query point is near to or within the bounding box of an internal node, recurse to its children.
    • If the query point is near to or within the bounding box of a leaf node, it's contribution to the GWN is computed exactly using a direct formula.

Implementation Details

  • Utilizing the existing BVH tree for this purpose required modifying it in the following ways:
    • Adding a LinearBVHTraverser::reduce_tree(<LeafAction>) method which evaluates the LeafAction lambda on each leaf node to produce a value, then adding those values up the tree so that each internal node is associated with the sum of values for all its children. This returns an array of these values for each internal node.
    • Adding an optional third positional argument to the search predicate lambda (i.e. the function which determines if a traverser should recurse at a given internal node) which contains an index to the internal node as laid out within the BVH.

This PR should consider the following before being un-drafted: Done!

  • Extend quest_winding_number_3d_ex to take in .stl files directly as input.
  • Add more command line options to the 2D and 3D examples to permit changing Taylor approximation order
  • Add more command line options to the 3D example to permit direct evaluation of a triangulated STEP file instead of fast approximation, even if this would be significantly slower without a difference in quality.
  • Add a few extra tests to make sure containment decisions between the current direct methods and fast approximation methods agree.
  • Add a method to the quest STEP reader to get a bounding box directly from the .STEP file using OpenCascade.

@jcs15c jcs15c self-assigned this Feb 2, 2026
@jcs15c jcs15c added enhancement New feature or request Quest Issues related to Axom's 'quest' component labels Feb 2, 2026
@jcs15c jcs15c marked this pull request as ready for review February 11, 2026 07:36
@jcs15c jcs15c changed the title Implement first fast approximation GWN methods Implement first set of fast approximation GWN methods Feb 11, 2026
public:
using NURBSCurve = axom::primal::NURBSCurve<double, 2>;
using CurveArrayView = axom::ArrayView<NURBSCurve>;
using CurveArrayView = axom::ArrayView<const NURBSCurve>;
Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I had to make this change for this to work more straight-forwardly with LinearizeCurves, since otherwise (to my knowledge) you can't use it to linearize a view of const objects without making a copy of all of them. But I think this is more close to how the class is actually implemented everywhere else, right? It doesn't seem to me like there's anywhere that a LinearizeCurves object ever needs to change a NURBSCurve in it's view.

@kennyweiss kennyweiss force-pushed the feature/kweiss/optimize-wn-3d branch from 9142cea to 4609f08 Compare February 14, 2026 01:32
@kennyweiss kennyweiss force-pushed the feature/kweiss/optimize-wn-3d branch from 4609f08 to cfd9c5c Compare February 25, 2026 20:36
Base automatically changed from feature/kweiss/optimize-wn-3d to develop March 3, 2026 16:46
@kennyweiss kennyweiss requested a review from publixsubfan March 4, 2026 19:18
@kennyweiss kennyweiss added this to the FY26 Summer release milestone Mar 4, 2026
@kennyweiss kennyweiss added Primal Issues related to Axom's 'primal component Spin Issues related to Axom's 'spin' component labels Mar 4, 2026
Copy link
Copy Markdown
Member

@kennyweiss kennyweiss left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Partial review: first batch of comments

Comment thread src/axom/primal/geometry/detail/analytic_test_surfaces.hpp Outdated
Comment thread src/axom/primal/geometry/detail/analytic_test_surfaces.hpp Outdated
Comment thread src/axom/primal/geometry/detail/analytic_test_surfaces.hpp Outdated
Comment thread src/axom/primal/geometry/detail/analytic_test_surfaces.hpp
Comment thread src/axom/primal/geometry/detail/analytic_test_surfaces.hpp
Comment thread src/axom/spin/policy/LinearBVH.hpp Outdated
}
else
{
// Do it in 2 stages.
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can you please expand on the two steps and why we need them for the non-seq policy?
Looks like step one is to compute on the leaf nodes and step 2 is to move it up the tree.
(whereas the seq policy branch does both)

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think the main benefit is just to improve the efficiency of operating over all the leaves, and I'm not certain that needs to be an exclusively parallel thing. I could probably rewrite reduce_recursion to not even take a lambda function, and just pass a view to the leaf array. I think @BradWhitlock originally wrote this, do you know of any particular reason it can't be done in two stages for both policies?

Comment thread src/axom/spin/internal/linear_bvh/bvh_traverse.hpp
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@jcs15c -- is there a simple example or unit test we can add to spin for the new reduce_tree traverser? E.g. adding a random collection of points to a BVH by their bounding boxes and the counting the number of points contained within each inner node?

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, I'm sure I can come up with something. One way that's agnostic to the way the tree is constructed is that, if you have 8 points added to the BVH, then 8 inner nodes should contain one point (the leaves), 4 contain two points, 2 contain four points, and 1 contains all eight points, right?

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Added!

Comment thread src/axom/primal/geometry/KnotVector.hpp Outdated
Copy link
Copy Markdown
Member

@kennyweiss kennyweiss left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks @jcs15c -- this looks great! I have a bunch of minor comments; nothing substantial.

Comment thread src/axom/quest/examples/quest_winding_number_2d.cpp Outdated
Comment thread src/axom/quest/examples/quest_winding_number_2d.cpp Outdated
Comment thread src/axom/quest/examples/quest_winding_number_3d.cpp Outdated
Comment thread src/axom/quest/io/STEPReader.hpp Outdated
std::string pjoin(const char *str, Args... args)
{
return axom::utilities::filesystem::joinPath(std::string(str), pjoin(args...));
}
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@BradWhitlock -- should we move these pjoin functions to core/utilities? Perhaps in a separate PR?

Comment thread src/axom/quest/FastApproximateGWN.hpp Outdated
* the collection of geometric objects.
*/
template <typename T, int NDIMS, int ORD>
class GWNMomentData
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Would it make sense to add some unit tests for the moment data computations? Or is this covered by the segment/triangle tests?

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Added some for 2D segments and 3D triangles!

Comment thread src/axom/quest/FastApproximateGWN.hpp Outdated
Comment thread src/axom/quest/FastApproximateGWN.hpp Outdated
Comment thread src/axom/quest/FastApproximateGWN.hpp
template <typename ExecSpace, typename ValueType, typename LeafAction>
axom::Array<ValueType> reduce_tree(LeafAction&& lf,
int allocatorID = axom::getDefaultAllocatorID()) const
{
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is the eventual plan for this function to be invokable from inside a GPU kernel as a device function? If not, I'd suggest moving it into the standard BVH interface.

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We'd like to, but are not there yet

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

As I understand it, a method like this couldn't go directly into the standard BVH interface because only LinearBVH describes how the data at internal nodes is actually arranged into arrays and arrayviews. This output of this function needs to align with that internal layout, and so if there was ever another implementation type for BVH, the actual logic of the function would have to change to reflect the internal layout of that other implementation. Does that seem consistent with how BVH.hpp and LinearBVH.hpp relate to one another?

@jcs15c jcs15c force-pushed the feature/spainhour/fast_gwn_triangles branch from 9993c59 to 0f31be3 Compare March 13, 2026 19:07
@jcs15c jcs15c force-pushed the feature/spainhour/fast_gwn_triangles branch from bc6e432 to 3f33b76 Compare March 16, 2026 17:15
Comment thread src/axom/quest/GWNMethods.hpp Outdated
poly_mesh,
AXOM_LAMBDA(axom::IndexType cellIdx,
const axom::numerics::Matrix<double>& coords,
[[maybe_unused]] const axom::IndexType* nodeIds) {
Copy link
Copy Markdown
Member

@kennyweiss kennyweiss Mar 16, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

hip is not happy w/ this attribute. Looks like we should change it to something like

Suggested change
[[maybe_unused]] const axom::IndexType* nodeIds) {
const axom::IndexType* /*nodeIds*/) {

... and similarly for the one below

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

No complaints here!

Copy link
Copy Markdown
Contributor

@Arlie-Capps Arlie-Capps left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@jcs15c , thank you for your hard work!

string(CONCAT _pass_regex
"INOUT_STATS:.*"
"pos_dofs=16.*"
"pos_dofs=14.*"
Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Just for future reference, the new OCC-derived AABB method has slightly different bounds, so there are different numbers of contained points.

@jcs15c jcs15c merged commit db9bead into develop Mar 23, 2026
15 checks passed
@jcs15c jcs15c deleted the feature/spainhour/fast_gwn_triangles branch March 23, 2026 19:44
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

enhancement New feature or request Primal Issues related to Axom's 'primal component Quest Issues related to Axom's 'quest' component Spin Issues related to Axom's 'spin' component

Projects

None yet

Development

Successfully merging this pull request may close these issues.

4 participants