Update 3D GWN method for trimmed NURBS surfaces #1566
Conversation
…ure/spainhour/winding_number_3d_update
…ure/spainhour/winding_number_3d_update
jcs15c
added
enhancement
| // Note: we want to find all intersections, so don't short-circuit | ||
| intersect_line_patch(line, | ||
| p1, | ||
| tp, | ||
| up, | ||
| vp, | ||
| order_u, | ||
| order_v, | ||
| u_offset, | ||
| u_scale, | ||
| v_offset, | ||
| v_scale, | ||
| sq_tol, | ||
| EPS, | ||
| isRay); | ||
|
|
||
| intersect_line_patch(line, | ||
| p2, | ||
| tp, | ||
| up, | ||
| vp, | ||
| order_u, | ||
| order_v, | ||
| u_offset + u_scale, | ||
| u_scale, | ||
| v_offset, | ||
| v_scale, | ||
| sq_tol, | ||
| EPS, | ||
| isRay); | ||
|
|
||
| intersect_line_patch(line, | ||
| p3, | ||
| tp, | ||
| up, | ||
| vp, | ||
| order_u, | ||
| order_v, | ||
| u_offset, | ||
| u_scale, | ||
| v_offset + v_scale, | ||
| v_scale, | ||
| sq_tol, | ||
| EPS, | ||
| isRay); | ||
|
|
||
| intersect_line_patch(line, | ||
| p4, | ||
| tp, | ||
| up, | ||
| vp, | ||
| order_u, | ||
| order_v, | ||
| u_offset + u_scale, | ||
| u_scale, | ||
| v_offset + v_scale, | ||
| v_scale, | ||
| sq_tol, | ||
| EPS, | ||
| isRay); | ||
| // *unless* we're already in a failure state | ||
| success = success && |
There was a problem hiding this comment.
This comment is a bit confusing since the calls are clearly being short-circuited.
Perhaps modify the comment to remind people that intersect_line_patch will return true unless the patch had degenerate intersections.
There was a problem hiding this comment.
It's getting to a length where I'd be tempted to use a loop over points, just to shrink the lines of code.
There was a problem hiding this comment.
I'm up for shortening this, especially since I think we'd want to go in the direction of adding an extra pass-by-reference argument to indicate a failure state instead of doing that with the return type.
EDIT: I've shortened up these calls, and in doing so, made the comments much more clear. The original idea was that we don't want to return early if intersections are found, which I don't think is super accurately described by short-circuiting anyway.
| // Lambda to generate a 3D rotation matrix from an angle and axis | ||
| // Formulation from https://en.wikipedia.org/wiki/Rotation_matrix#Axis_and_angle | ||
| auto angleAxisRotMatrix = [](double theta, const Vector<T, 3>& axis) -> numerics::Matrix<T> { | ||
| const auto unitized = axis.unitVector(); | ||
| const double x = unitized[0], y = unitized[1], z = unitized[2]; | ||
| const double c = cos(theta), s = sin(theta), C = 1 - c; | ||
|
|
||
| auto matx = numerics::Matrix<T>::zeros(3, 3); | ||
|
|
||
| matx(0, 0) = x * x * C + c; | ||
| matx(0, 1) = x * y * C - z * s; | ||
| matx(0, 2) = x * z * C + y * s; | ||
|
|
||
| matx(1, 0) = y * x * C + z * s; | ||
| matx(1, 1) = y * y * C + c; | ||
| matx(1, 2) = y * z * C - x * s; | ||
|
|
||
| matx(2, 0) = z * x * C - y * s; | ||
| matx(2, 1) = z * y * C + x * s; | ||
| matx(2, 2) = z * z * C + c; | ||
|
|
||
| return matx; | ||
| }; |
There was a problem hiding this comment.
In case it's helpful, @BradWhitlock recently added several utility functions for generating transformations: https://github.com/LLNL/axom/blob/develop/src/axom/core/numerics/transforms.hpp
Perhaps we should move this angle-axis there? It appears a few times in axom
(No need to make this change for the current PR)
There was a problem hiding this comment.
I think that makes sense, especially since the interface is very similar to the ones in that file. Not a huge change either way, and keeps this already pretty meaty function more consice.
| // ...unless the average direction is zero | ||
| double theta = axom::utilities::random_real(0.0, 2 * M_PI); | ||
| double u = axom::utilities::random_real(-1.0, 1.0); | ||
| cast_direction = Vector<T, 3> {sin(theta) * sqrt(1 - u * u), cos(theta) * sqrt(1 - u * u), u}; |
There was a problem hiding this comment.
Is it worth computing sqrt(1 - u * u) one time here?
There was a problem hiding this comment.
It probably is technically faster since there's a call to sqrt involved, but this branch of this if-stament is the far less common one, and the total cost is completely dwarfed by the rest of the algorithm.
…om:LLNL/axom into feature/spainhour/winding_number_3d_update
…om:LLNL/axom into feature/spainhour/winding_number_3d_update
Summary
Altogether, this significantly increases the processing speed of points which are close to the surface. Points which are close to a boundary of the input surface are evaluated more accurately with the new method, but at a comparable computational expense. Points which are far from the boundary are treated identically to the previous version of the algorithm.
This addition is facilitated by
NURBSPatch::calculateUntrimmedPatchNormalmethod which defines the direction of the line used in the intersection routine. This utilizes newevaluate_integral(NURBSCurve)methods, all of which come with associated tests.intersect(NURBSPatch/BezierPatch, Line/Ray)routines so that the return value indicates the success of the algorithm (i.e., the algorithm does not return early due to degenerate intersections). Previously, this return value indicated if intersections were found, which is equivalently found by checking the size of the return vectors.Additional Changes
primal::intersect_circle_bboxso that relevant routines can more robustly handle floating point errors, addressing issues Re-enableprimal_intersecttest #1572 and Add tolerance to primal's intersect_circle_bbox #1573.