Add support for objects specified by NURBS?

[elbert5770]

I've worked through the dogfish example, which produces an sdf from a B-spline and is very nice. However, it's my understanding that the sdf generating function doesn't measure closest distance but are just conformal maps in the y-direction? Will this cause some numerical error when calculating the gradients? The projection (and closest distance) of any point onto the curve should be orthogonal to the curve.

If the order and knot vector of the NURBS or B-spline is such that second derivatives exist, then equation 6.3 in the NURBS Book 2nd ed by Piegl and Tiller allows calculation of the sdf. The downside is that you need to check every span. Additionally, if the curve is not closed or C2 continuous, it will make everything so much more complicated. That said, the continuity can be simply read off the knot vector. Additionally, it makes physical sense that the head and tail of the dogfish are C0 continuous.

Anyway, I think this more of a topic for discussion rather than an action item.

I've searched for Python, Julia and C++ tools to calculate the distance of a point from a NURBS curve and haven't found anything on github yet. Of the NURBS packages in Julia, the BasicBSpline package is in active development and is pretty well organized. It's just one person though and while off to a good start, there's still a lot left to do. Python-NURBS is really good for curve fitting.

[weymouth]

Good idea Don. I agree that this should start with a discussion.

You are correct that the method in the dogfish shark example doesn't generate a true SDF to the shark. It warps a line segment SDF in two ways that don't maintain the minimum distance property.

It's worth noting that we don't need a perfect SDF for WaterLily. The values a small distance from the fluid/body interface are clamped

https://github.com/weymouth/WaterLily.jl/blob/af00631fa7b1fe1099c87360e5f39df4b2dd7724/src/Body.jl#L26-L27

and we scale by the AD gradient to improve the distance function within those tight bounds

https://github.com/weymouth/WaterLily.jl/blob/af00631fa7b1fe1099c87360e5f39df4b2dd7724/src/AutoBody.jl#L96-L98

However, your point stands. The direction of the normal and the velocity are not corrected, and large amplitude warping will give very bad values and unphysical simulations.

[weymouth]

Also, functions that compute the SDF (and normal and velocity) to a NURBS should really live in an external package. They aren't WaterLily specific at all. Do you think the https://github.com/hyrodium/BasicBSpline.jl package is the place to do this?

I don't have any NURBS code, but I have previously worked out a hacky method to get the SDF for a bspline:

using Roots
function bezier_sdf(P0,P1,P2,P3)
    # Define cubic Bezier
    A,B,C,D = P0,-3P0+3P1,3P0-6P1+3P2,-P0+3P1-3P2+P3
    curve(t) = A+B*t+C*t^2+D*t^3
    dcurve(t) = B+2C*t+3D*t^2

    # Signed distance to the curve
    candidates(X) = union(0,1,find_zeros(t -> (X-curve(t))'*dcurve(t),0,1,naive=true,no_pts=3,xatol=0.01))
    function distance(X,t)
        V = X-curve(t)
        D = dcurve(t)
        copysign(√(V'*V),V[2]*D[1]-V[1]*D[2])
    end
    X -> argmin(abs, distance(X,t) for t in candidates(X))
end

This isn't optimized at all, and I would love to see a better approach that works on more general splines.

[weymouth]

I generalized my code above, but needed to change it a lot to run on GPUs since Roots.jl doesn't seems to work on CUDA.

Here's the updated version:

using StaticArrays,ForwardDiff
"""
    parametric_sdf(curve)

Create a function to find the signed distance from a point `X` in 2D space
to a parametric `curve(t)` where `t ∈ [0,1]`. The resulting function uses 
ForwardDiff and a simple multi-start Newton root-finding method to identifiy 
`t` closest to `X`. The returned function runs on GPU and has optional arguments:

    function(X,tol=0.01,t₀=(0.01,0.5,0.99),itmx=5)

where `tol` is the tolerance in `t`, `t₀` are the multi-start points (along 
with the bounds 0,1), and `itmx` is the max iterations. For example:

    R = 2
    circle_SDF = parametric_sdf(t->SA[R*sin(2π*t),R*cos(2π*t)])
    @assert circle_SDF(SA[3,4],tol=1e-8) ≈ 5-R

Note that sign of the SDF depends on the direction of `curve`. The example above
wraps clockwise such that distances outside the circle are positive.
"""
function parametric_sdf(curve)
    dcurve(t) = ForwardDiff.derivative(curve,t)
    ddcurve(t) = ForwardDiff.derivative(dcurve,t)
    norm(t) = (s=dcurve(t); SA[-s[2],s[1]])
    vector(X,t) = X-curve(t)
    align(X,t) = vector(X,t)'*dcurve(t)
    dalign(X,t) = vector(X,t)'*ddcurve(t)-sum(abs2,dcurve(t))
    distance(X,t) = (v=vector(X,t); copysign(√(v'*v),norm(t)'*v))
    function(X; tol=0.01,t₀=(0.01,0.5,0.99),itmx=5)
        # distance to ends
        dmin = distance(X,0.)
        d = distance(X,1.)
        abs(d)<abs(dmin) && (dmin=d)
        # check for smaller distance along curve
        for t in t₀
            for _ in 1:itmx # Newton-method to find where align(t)=0 
                dt = align(X,t)/dalign(X,t)
                t -= dt
                t = clamp(t,0.,1.)
                (abs(dt)<tol || t==0 || t==1) && break
            end
            d = distance(X,t)
            abs(d)<abs(dmin) && (dmin=d)
        end
        dmin
    end
end

[elbert5770]

This means Waterlily.jl can handle any arbitrary shape as long as you define the function to traverse the curve clockwise? Very cool. I will check this out.

[weymouth]

This example isn't in the code yet. And this is only for 1D curves, and there is no NURBS support.

[weymouth]

I've created ParametricBodies.jl as a separate repo which further generalizes this code. It's still only for 1D curves, but further development discussion should take place over there.