Accelerates 2D/3D GWN queries with memoization

[jcs15c]

This PR is ready to go (assuming cooperation with the Gitlab CI), but I'm waiting to merge it until I'm back in office on 9/30, since it isn't critical for the upcoming release.

Summary

The current 2D and 3D algorithms currently compute redundant information between evaluation at different query points, namely evaluations of the curve/surface. This PR is an algorithm refactoring/enhancement that addresses this by utilizing memoization for performance improvements. A NURBSCurveGWNCache or NURBSPatchGWNCache object stores the original NURBS object along with map(s) that stores these intermediate values as they are used in the GWN algorithm:

• In 2D, the input NURBS curve is decomposed into Bezier curves, and necessary levels of subdivisions are stored for each Bezier curve. This implementation is considerably simpler than before, e.g. no longer constructs and checks polygons at each refinement, but avoids extra memory allocations.
• In 3D, quadrature nodes and tangent vectors to the surface are stored for each trimming curve in a patch. Furthermore, a copy of the original patch altered to improve the robustness of the GWN algorithm, e.g. linearly extending the parameter space of the patch so that the untrimmed geometry is larger than the visible region.

These methods are compatible with user supplied cache objects, but are also constructed automatically in cases where the GWN algorithm is called on arrays of query points to simplify use.

Note that these acceleration techniques are currently not reflected in the winding number sampling shaper code in quest, namely because these methods are built on CurvedPolygon objects which contain Bezier, not NURBS curves. This will be addressed in a future PR.

Some lingering questions that can hopefully be thought about during review:

• Is this framework of dedicated "cache objects" the implementation that is most conducive to future changes, e.g., GPU implmentation?
• Is GWNCache the correct label for these objects? Conceivably, there could be other kinds of memoization strategies that involve curves/surfaces, but other geometric operations, such as closest point.

Performance Improvements

To evaluate the performance improvement of the 3D memoization scheme, I've executed this block of code on this branch and develop on my local laptop. On a release build, this code executed in 13.8 minutes on develop and 10.3 seconds on this branch.

TEST(primal_integral, demo_test)
{
  using Point3D = primal::Point<double, 3>;
  using Vector3D = primal::Vector<double, 3>;

  // Set up an array of points to test against
  axom::Array<Point3D> points;

  for(int i = 0; i < 1e5; ++i)
  {
    auto x = axom::utilities::random_real( -2.0, 2.0 );
    auto y = axom::utilities::random_real( -2.0, 2.0 );
    auto z = axom::utilities::random_real( -2.0, 2.0 );
    points.push_back( Point3D { x, y, z } );
  }
  
  // Get the array of points
  auto sphere_faces = make_sphere_biquintic();

  // Evaluate the winding number at the points
  auto gwn = winding_number(points, sphere_faces);
}

For 2D, we run this block of code, which takes 108ms on the current develop branch and 18ms on the current branch.

TEST(primal_integral, demo_test)
{
  primal::Point<double, 2> data[7] = {primal::Point<double, 2> {1.0, 0.0},
                                      primal::Point<double, 2> {1.0, 2.0},
                                      primal::Point<double, 2> {-1.0, 2.0},
                                      primal::Point<double, 2> {-1.0, 0.0},
                                      primal::Point<double, 2> {-1.0, -2.0},
                                      primal::Point<double, 2> {1.0, -2.0},
                                      primal::Point<double, 2> {1.0, 0.0}};
  double weights[7] = {1.0, 1. / 3., 1. / 3., 1.0, 1. / 3., 1. / 3., 1.0};
  double knots[11] = {0.0, 0.0, 0.0, 0.0, 0.5, 0.5, 0.5, 1.0, 1.0, 1.0, 1.0};

  primal::NURBSCurve<double, 2> circle(data, weights, 7, knots, 11);
  //primal::detail::NURBSCurveGWNCache circle_data(circle);

  // Set up an array of points to test against
  axom::Array<primal::Point<double, 2>> points;

  for(int i = 0; i < 1e5; ++i)
  {
    auto x = axom::utilities::random_real(-2.0, 2.0);
    auto y = axom::utilities::random_real(-2.0, 2.0);
    points.push_back(primal::Point<double, 2> {x, y});
  }

  // Evaluate the winding number at the points
  for(auto & pt : points)
    winding_number(pt, circle_data); // `circle` on develop
}

Additional changes and bugfixes

• Updates the quest_winding_number.cpp example to utilize cached NURBS curves.
• An additional NURBSPatch::expandParameterSpace method that allows the patch to the expanded differently along each axis.
• Changes several internal intersection tolerances for consistency.
• Adds isHalfOpen option for ray-bezier intersection tests.
• Adds a method to evaluate average surface normals for trimmed surfaces.

[kennyweiss]

Really nice improvements @jcs15c and great performance numbers!

Please update the RELEASE_NOTES.

[BradWhitlock]

Nice speedups.

[Arlie-Capps]

isOnCurve should be [out] or [inout]

[Arlie-Capps]

I suspect all the Doxygen should be reviewed for this.

[jcs15c]

Very true. There's a couple of places where the [in]/[out]/[inout] labels are missing altogether. I'll go through and double-check everything.