Adds untrimmed NURBS patches to primal

[jcs15c]

Summary

• This PR is a feature
• It adds a NURBSPatch class to primal. NURBSPatch has the following features, along with relevant tests:
  • Construictors from arrays of points (optionally weights) and either the degree or knot vector
  • Evaluation/Derivative methods
  • Knot insertion
  • Patch subdivision/Bezier extraction
• It also enhances the Bezier extraction method of NURBSCurve to be more efficient, performing the knot insertion in place to avoid repeated calculations.
  • An extra test has been added to primal_nurbs_curve.cpp to test a currently unchecked edge case.
• Minor bug-fixes in NURBSCurve and KnotVector

Additional Changes

• Rework the KnotVector and NURBSCurve classes to be more tolerant of parameters slightly outside the range of the knots.

[kennyweiss]

Thank you @jcs15c for this great new feature w/ comprehensive tests!

[kennyweiss]

Are there cases where the derivative is not defined, e.g. near cusps? If so, what will the normal be?

[jcs15c]

It's possible that one or both derivatives are zero, in which case the normal would be zero too. That doesn't mean that the surface itself doesn't have a normal there (like for a corner of a degenerate triangular patch), but the normal with respect to this parameterization is zero. Accounting for all the relevant edge cases is tricky, and I haven't worked out a consistent system for that yet.

[kennyweiss]

Would it make sense to add a note? Or perhaps later when the edge cases are more worked out?

[jcs15c]

Couldn't hurt to add a comment here, although it isn't like returning a zero vector is itself incorrect. For a simple cusp on a 2D curve, the tangent vector at the point of the cusp really is just continuously zero, and although I have less intuition for the 3D case, I'm fairly sure that's the case there too. Even though the direction of the tangent/normal vector changes depending on the direction of approach, the magnitude goes to zero regardless. In my experience, the more pressing issue is when you need to take a unit normal vector at a point like that, since that genuinely is always going to be discontinuous at the cusp. This is a bit troublesome with the current primal::Vector implementation, which returns something like (1,0,0) by default if the norm is zero.

As an example of how I've dealt with these geometric/unit vectors before, the 2D winding number classes need to know tangent vectors at endpoints, and in those cases I manually iterated over control points until I found a pair that correctly defined a geometric tangent vector. That doesn't strictly help with the cusp case, but in 2D you could imagine splitting the curve at the point of discontinuity and handling each half that way.

[kennyweiss]

Thanks for the updates!

[bmhan12]

LGTM!

One thing missing is NURBSPatch.hpp should be added to the list of headers in primal/CMakeLists.txt

[jcs15c]

I made a few extra changes to the KnotVector class and the NURBS objects to make them more tolerant of input parameter values that are slightly outside of the range of the knot values, which happens surprisingly often with floating point arithmetic. These checks for "slightly outside the range" are handled by the KnotVector::isValidParameter method, which is the only place where the exact value of this numerical tolerance is kept.

In the future, it might be desirable to give NURBS curves and surfaces the ability to extrapolate beyond their knot spans slightly, as is naturally the case for the analogous Bezier objects, but that's a comparatively computationally expensive procedure compared to simple clamping to handle floating point errors. This is also why I don't think it's necessary to expose this tolerance parameter to a user, and to some extent, I prefer the idea of all instances of these objects using the same tolerance level.

@kennyweiss What do you think of this as a solution?

[kennyweiss]

I made a few extra changes to the KnotVector class and the NURBS objects to make them more tolerant of input parameter

@kennyweiss What do you think of this as a solution?

Thanks @jcs15c -- I think this is a reasonable and sufficient solution for now. If we need to extrapolate at parametric points near the boundary in the future, it should be easy to add that capability (at the expense of some additional computational work, as you mentioned).