Issues with SDF of complex dynamic NURBS body

[evanvanderweide]

The Issue

I am encountering weird behaviour in the signed distance field (SDF) generated from a DynamicNurbsBody in ParametricBodies.jl. During simulations using a jellyfish-shaped body, random geometry cells appear inside the wake field, despite no geometry being present there. Additionally, near regions with sharp corners, the SDF includes an extra geometry cell.

Any guidance on whether this is a bug in the SDF implementation or something incorrect in my implementation would be appreciated.

Results
I've included some visuals that illustrate this problem below. The first figure shows the values of the SDF generated using the minimal working example provided below, with some weird geometry cells present. The second one is a time frame from a full simulation of the jellyfish with a sim.L = 2^7. It shows how the generated cells interact with my wake field. As a result, this impacts the results that I compute with the simulation.

For a Nurbs Curve with a degree of 3, the problem becomes even worse:

What I Tried
I varied several parameters:

• grid size
• number of control points,
• adjusting ϵ
• the polynomial degree of the Nurbs Curve
• assigning a thickness to the body

In all cases, the problems kept occurring; only the extent of the problems would change. For a polynomial with degree 2 as a curve, the problems seem minimal, but there is still a dot at the corners of the geometry and a couple of ghost cells in the wake.

Minimal Working Example
See below a minimal working example that can be used to generate the SDFs of the jellyfish geometry. I define four control point sets here that are problematic, picked from the full jellyfish motion cycle. This code uses the SDF function defined in ParametricBodies, but this yields the same results as the WaterLily function.

I am using the following package versions:

• WaterLily: v1.5.2
• ParametricBodies: v0.3.2

using LinearAlgebra, Plots, ParametricBodies, StaticArrays, WaterLily
knots_vector(p::Int, Ncp::Int) = vcat(zeros(p+1), (Ncp-p-1 > 0 ? collect(range(0.0, 1.0, length=Ncp-p+1))[2:end-1] : Float64[]), ones(p+1))

ϵ = 1; deg = 3; D = 2^7; xs = range(0, 6D, step=1); ys = range(0, 6D, step=1) 

cps_0 = SA{Float64}[256.0 256.0 256.0 256.0 257.263 258.619 260.095 261.738 263.626 265.942 269.131 273.131 277.703 282.723 288.121 293.934 300.147 306.724 313.647 320.882 328.413 336.235 344.288 352.458 360.677 368.911 377.134 385.324 393.463 401.532 409.524 417.444 425.295 433.094 440.891 448.711 456.566 464.471 472.442 480.457 488.334 495.37 499.374 501.598 503.29 504.727 505.322 505.937 501.368 499.682 498.111 496.03 491.565 483.737 475.546 467.338 459.133 450.912 442.699 434.488 426.279 418.064 409.849 401.637 393.455 385.342 377.408 369.819 362.658 355.938 349.718 344.003 338.9 335.113 332.341 330.014 327.944 326.045 327.944 330.014 332.341 335.113 338.9 344.003 349.718 355.938 362.658 369.819 377.408 385.342 393.455 401.637 409.849 418.064 426.279 434.488 442.699 450.912 459.133 467.338 475.546 483.737 491.565 496.03 498.111 499.682 501.368 505.937 505.322 504.727 503.29 501.598 499.374 495.37 488.334 480.457 472.442 464.471 456.566 448.711 440.891 433.094 425.295 417.444 409.524 401.532 393.463 385.324 377.134 368.911 360.677 352.458 344.288 336.235 328.413 320.882 313.647 306.724 300.147 293.934 288.121 282.723 277.703 273.131 269.131 265.942 263.626 261.738 260.095 258.619 257.263 256.0 256.0 256.0 256.0; 320.0 317.44 314.88 312.32 311.855 303.732 295.629 287.556 279.542 271.644 264.068 256.872 250.025 243.497 237.277 231.446 226.04 221.082 216.619 212.691 209.388 206.879 205.213 204.198 203.718 203.7 204.115 204.94 206.157 207.76 209.69 211.88 214.286 216.858 219.444 221.95 224.279 226.134 227.73 229.134 231.588 235.304 242.206 250.075 258.038 266.022 273.798 282.059 280.509 273.054 265.192 257.377 250.789 249.028 248.481 248.185 248.029 247.748 247.496 247.288 247.139 247.074 247.135 247.397 248.047 249.346 251.478 254.624 258.656 263.37 268.738 274.657 281.102 288.372 296.119 304.012 311.982 320.0 328.018 335.988 343.881 351.628 358.898 365.343 371.262 376.63 381.344 385.376 388.522 390.654 391.953 392.603 392.865 392.926 392.861 392.712 392.504 392.252 391.971 391.815 391.519 390.972 389.211 382.623 374.808 366.946 359.491 357.941 366.202 373.978 381.962 389.925 397.794 404.696 408.412 410.866 412.27 413.866 415.721 418.05 420.556 423.142 425.714 428.12 430.31 432.24 433.843 435.06 435.885 436.3 436.282 435.802 434.787 433.121 430.612 427.309 423.381 418.918 413.96 408.554 402.723 396.503 389.975 383.128 375.932 368.356 360.458 352.444 344.371 336.268 328.145 327.68 325.12 322.56 320.0]
cps_1 = SA{Float64}[256.0 256.0 256.0 256.0 257.255 258.601 260.066 261.693 263.56 265.843 268.987 272.951 277.494 282.492 287.869 293.667 299.867 306.434 313.35 320.575 328.094 335.907 343.957 352.127 360.35 368.59 376.819 385.019 393.171 401.257 409.27 417.214 425.094 432.921 440.746 448.594 456.475 464.403 472.397 480.431 488.353 495.359 499.019 500.849 502.132 503.155 503.198 503.591 499.299 498.145 496.983 495.274 490.997 483.132 474.93 466.712 458.497 450.267 442.046 433.827 425.612 417.393 409.179 400.973 392.801 384.705 376.792 369.23 362.105 355.436 349.263 343.584 338.518 334.808 332.104 329.839 327.826 325.984 327.826 329.839 332.104 334.808 338.518 343.584 349.263 355.436 362.105 369.23 376.792 384.705 392.801 400.973 409.179 417.393 425.612 433.827 442.046 450.267 458.497 466.712 474.93 483.132 490.997 495.274 496.983 498.145 499.299 503.591 503.198 503.155 502.132 500.849 499.019 495.359 488.353 480.431 472.397 464.403 456.475 448.594 440.746 432.921 425.094 417.214 409.27 401.257 393.171 385.019 376.819 368.59 360.35 352.127 343.957 335.907 328.094 320.575 313.35 306.434 299.867 293.667 287.869 282.492 277.494 272.951 268.987 265.843 263.56 261.693 260.066 258.601 257.255 256.0 256.0 256.0 256.0; 320.0 317.44 314.88 312.32 311.847 303.718 295.607 287.525 279.5 271.586 263.986 256.765 249.893 243.342 237.098 231.242 225.813 220.833 216.347 212.387 209.042 206.483 204.778 203.726 203.208 203.152 203.524 204.299 205.46 206.998 208.857 210.974 213.305 215.804 218.322 220.762 223.009 224.715 226.164 227.437 229.804 233.48 240.563 248.628 256.767 264.913 272.779 281.197 279.434 271.817 263.749 255.723 249.01 247.406 246.949 246.751 246.701 246.481 246.292 246.154 246.081 246.102 246.26 246.635 247.393 248.789 251.005 254.222 258.324 263.112 268.539 274.497 280.974 288.286 296.06 303.975 311.964 320.0 328.036 336.025 343.94 351.714 359.026 365.503 371.461 376.888 381.676 385.778 388.995 391.211 392.607 393.365 393.74 393.898 393.919 393.846 393.708 393.519 393.299 393.249 393.051 392.594 390.99 384.277 376.251 368.183 360.566 358.803 367.221 375.087 383.233 391.372 399.437 406.52 410.196 412.563 413.836 415.285 416.991 419.238 421.678 424.196 426.695 429.026 431.143 433.002 434.54 435.701 436.476 436.848 436.792 436.274 435.222 433.517 430.958 427.613 423.653 419.167 414.187 408.758 402.902 396.658 390.107 383.235 376.014 368.414 360.5 352.475 344.393 336.282 328.153 327.68 325.12 322.56 320.0]
cps_2 = SA{Float64}[256.0 256.0 256.0 256.0 257.246 258.58 260.03 261.639 263.483 265.732 268.827 272.751 277.263 282.236 287.592 293.375 299.566 306.129 313.042 320.26 327.77 335.576 343.626 351.801 360.03 368.279 376.521 384.734 392.904 401.012 409.051 417.025 424.937 432.8 440.657 448.538 456.449 464.405 472.424 480.477 488.436 495.396 498.684 500.108 500.973 501.574 501.094 501.194 497.184 496.532 495.778 494.438 490.346 482.44 474.223 465.992 457.763 449.52 441.287 433.056 424.831 416.605 408.388 400.187 392.026 383.946 376.054 368.519 361.432 354.816 348.697 343.059 338.037 334.415 331.786 329.589 327.641 325.861 327.641 329.589 331.786 334.415 338.037 343.059 348.697 354.816 361.432 368.519 376.054 383.946 392.026 400.187 408.388 416.605 424.831 433.056 441.287 449.52 457.763 465.992 474.223 482.44 490.346 494.438 495.778 496.532 497.184 501.194 501.094 501.574 500.973 500.108 498.684 495.396 488.436 480.477 472.424 464.405 456.449 448.538 440.657 432.8 424.937 417.025 409.051 401.012 392.904 384.734 376.521 368.279 360.03 351.801 343.626 335.576 327.77 320.26 313.042 306.129 299.566 293.375 287.592 282.236 277.263 272.751 268.827 265.732 263.483 261.639 260.03 258.58 257.246 256.0 256.0 256.0 256.0; 320.0 317.44 314.88 312.32 311.836 303.694 295.571 287.476 279.436 271.502 263.872 256.62 249.717 243.137 236.862 230.979 225.524 220.522 216.013 212.019 208.629 206.017 204.266 203.174 202.616 202.52 202.847 203.569 204.667 206.133 207.915 209.953 212.203 214.624 217.07 219.441 221.61 223.194 224.528 225.712 228.016 231.696 238.959 247.195 255.483 263.766 271.712 280.225 278.172 270.403 262.156 253.945 247.114 245.641 245.256 245.14 245.178 245.017 244.893 244.826 244.836 244.952 245.224 245.73 246.616 248.126 250.438 253.737 257.917 262.791 268.29 274.296 280.816 288.179 295.986 303.928 311.942 320.0 328.058 336.072 344.014 351.821 359.184 365.704 371.71 377.209 382.083 386.263 389.562 391.874 393.384 394.27 394.776 395.048 395.164 395.174 395.107 394.983 394.822 394.86 394.744 394.359 392.886 386.055 377.844 369.597 361.828 359.775 368.288 376.234 384.517 392.805 401.041 408.304 411.984 414.288 415.472 416.806 418.39 420.559 422.93 425.376 427.797 430.047 432.085 433.867 435.333 436.431 437.153 437.48 437.384 436.826 435.734 433.983 431.371 427.981 423.987 419.478 414.476 409.021 403.138 396.863 390.283 383.38 376.128 368.498 360.564 352.524 344.429 336.306 328.164 327.68 325.12 322.56 320.0]
cps_3 = SA{Float64}[256.0 256.0 256.0 256.0 257.234 258.554 259.988 261.577 263.396 265.609 268.653 272.534 277.012 281.959 287.294 293.063 299.248 305.81 312.725 319.94 327.446 335.248 343.301 351.483 359.723 367.985 376.242 384.474 392.665 400.8 408.87 416.879 424.83 432.732 440.629 448.546 456.493 464.48 472.525 480.598 488.585 495.481 498.374 499.383 499.826 500.001 499.04 498.777 495.049 494.859 494.507 493.529 489.616 481.667 473.432 465.183 456.937 448.676 440.428 432.182 423.942 415.705 407.483 399.285 391.133 383.07 375.199 367.69 360.641 354.081 348.022 342.43 337.46 333.935 331.388 329.265 327.388 325.675 327.388 329.265 331.388 333.935 337.46 342.43 348.022 354.081 360.641 367.69 375.199 383.07 391.133 399.285 407.483 415.705 423.942 432.182 440.428 448.676 456.937 465.183 473.432 481.667 489.616 493.529 494.507 494.859 495.049 498.777 499.04 500.001 499.826 499.383 498.374 495.481 488.585 480.598 472.525 464.48 456.493 448.546 440.629 432.732 424.83 416.879 408.87 400.8 392.665 384.474 376.242 367.985 359.723 351.483 343.301 335.248 327.446 319.94 312.725 305.81 299.248 293.063 287.294 281.959 277.012 272.534 268.653 265.609 263.396 261.577 259.988 258.554 257.234 256.0 256.0 256.0 256.0; 320.0 317.44 314.88 312.32 311.82 303.663 295.523 287.41 279.35 271.392 263.728 256.438 249.498 242.881 236.571 230.656 225.174 220.148 215.617 211.589 208.152 205.484 203.682 202.544 201.944 201.808 202.089 202.753 203.783 205.17 206.868 208.823 210.989 213.329 215.699 217.998 220.095 221.59 222.849 223.991 226.26 229.99 237.43 245.803 254.207 262.594 270.607 279.154 276.735 268.825 260.435 252.072 245.132 243.758 243.423 243.368 243.474 243.368 243.308 243.315 243.413 243.632 244.032 244.689 245.723 247.362 249.782 253.171 257.438 262.407 267.989 274.055 280.628 288.05 295.896 303.871 311.914 320.0 328.086 336.129 344.104 351.95 359.372 365.945 372.011 377.593 382.562 386.829 390.218 392.638 394.277 395.311 395.968 396.368 396.587 396.685 396.692 396.632 396.526 396.632 396.577 396.242 394.868 387.928 379.565 371.175 363.265 360.846 369.393 377.406 385.793 394.197 402.57 410.01 413.74 416.009 417.151 418.41 419.905 422.002 424.301 426.671 429.011 431.177 433.132 434.83 436.217 437.247 437.911 438.192 438.056 437.456 436.318 434.516 431.848 428.411 424.383 419.852 414.826 409.344 403.429 397.119 390.502 383.562 376.272 368.608 360.65 352.59 344.477 336.337 328.18 327.68 325.12 322.56 320.0]
cps_set = [cps_0, cps_1, cps_2, cps_3]

plt = Plots.plot(; aspect_ratio=:equal, legend=:topright,title="Jellyfish geometry evolution", xlabel="x", ylabel="y")
for (i,cps) in enumerate(cps_set)
    body        =   DynamicNurbsBody(NurbsCurve(cps, (knots_vector(deg, size(cps, 2))), ones(Float64, size(cps, 2)));thk=0, boundary=true)
    vals        =   [sdf(body, WaterLily.loc(0, CartesianIndex(x, y), Float64)) for y in ys, x in xs]
    sdf_plot    =   Plots.heatmap(xs, ys, vals; color=:viridis, xlims=(0,6D),ylims=(0,6D), aspect_ratio=1, title="Signed Distance Field $i for degree $deg")
                    Plots.contour!(xs, ys, vals, linewidth=2, color=:leonardo, levels=[0])  # Contour where sdf=0
    display(sdf_plot)
end

Thanks a lot in advance for the response.

[weymouth]

This is certainly possible. My test cases for the NURBS curve locator don't include anything with extremely tight rounded corners like you have here.

[weymouth]

Yeah, there are two things going on. First off, you have way more control points than you need for a degree=3 curve. By far the easiest fix is to use degree=1 (a linear spline) for these shapes. This gives you a fast/correct SDF. LilyPad used linearly splines for everything, so I don't think there will be any issue.

for (i,cps) in enumerate(cps_set)
    body     = ParametricBody(BSplineCurve(cps,degree=deg))
    vals     = [sdf(body, SA[x, y], 0.) for y in ys, x in xs]
    sdf_plot = Plots.heatmap(xs, ys, vals; color=:viridis, xlims=(0,6D),ylims=(0,6D), aspect_ratio=1, title="Signed Distance Field $i for degree $deg")
    Plots.contour!(xs, ys, vals, linewidth=2, color=:leonardo, levels=[0])  # Contour where sdf=0
    plot!(body.curve)
    display(sdf_plot)
end

It also looks like there might be an issue with your control points near the very top of the bell (left side of the body in the image aboe). There seem to be a large number of overlapping points there. You might get better results if you clean that up.

Note that there is a problem/bug with the NurbsLocator given an extremely non-convex shape like this and degree>1. The degree=1 minimizer set the initial condition for the iterative locator and that guess can be extremely bad for shapes like this. I'll think about how to fix this...

[evanvanderweide]

Thank You for the response.

I tried it with ParametricBody(BSplineCurve(cps,degree=1)) and a reduced number of control points, and it does indeed help when generating the SDF plots.

For some other frames in the full jellyfish motion cycle, some issues occurred in the SDF, but I managed to minimise these by adjusting the angles between control points at the flap edge. It now behaves better this way compared to my previous versions. I will simulate with this on a fine grid to capture what happens.

EDIT: I ran a simulation of the full motion cycle, and it behaves well, so this is a solution I can continue with for now.

[weymouth]

Nice! The time is O(N_control * N_grid), so reducing the control points will be faster. I can't think of a reason that reducing points would make the distance function more "stable" for a linear spline. If you have an example of that, please post it, so I can check.

Also, you should be able to use DynamicNurbsBody, just with a degree=1 curve. I was just cleaning up this example since the problem you reported was with the SDF, not the velocity. Have you check the velocity field?

[evanvanderweide]

I agree with you about the control points. When I change the number on the control points with my new body definition, it indeed does not change. So the fact that it works now results from turning it into a degree=1 curve and my adjustments to the flap definition. And indeed, when I change it to a degree=2 or degree=3 curve again, problems are present, so that confirms it.

Something I figured out this morning (because You said I should be able to use the DynamicNurbsBody I tried to use this again) was that having my knots vector as Float32 types messes up the SDF in my case. I noticed this as DynamicNurbsBody(NurbsCurve(cps,knots,weights)) did not work and ParametricBody(BSplineCurve(cps; degree=deg)) did work. NurbsCurve(cps,knots,weights) ended up working with a Float64 type knots vector. I plotted both versions below (sorry for the basic plot, but it's simply the contour of SD = 0). I have varied the Type of my control points but never for my knots, so it could have benefited my SDF problems.

The code that I used to produce this:

knots_vector(p::Int, Ncp::Int) = vcat(zeros(p+1), (Ncp-p-1 > 0 ? collect(range(0.0, 1.0, length=Ncp-p+1))[2:end-1] : Float64[]), ones(p+1))
@inline function TwoDimJellyfish(::Type{T}=Float32; D=2^7, Re=302, U=1, ϵ=1, thk=0, deg=1, mem=Array) where {T<:AbstractFloat}
    ν           =   U * D / Re

    cps         =   SA[256.0 256.0 256.0 256.0 257.62757550268054 259.2551510053611 261.06528451412305 264.78561568349545 268.50848794620987 276.3700768839923 284.54018037667265 293.37207429232427 304.3929798456619 315.4237870831155 327.433951108558 340.2637498776179 353.1494974328214 366.4794486555612 379.8213764417027 393.13929540174865 406.44017404365934 419.65751863205605 432.69379302804185 445.72953530402515 458.7464922215454 471.75749732833503 484.76366010020786 496.3233469543589 500.2611876427905 500.3816528542693 500.4547776492527 500.5111761579894 493.01617917024464 491.8619904979994 491.4209636271727 490.79624246528175 481.8841018259439 468.55207673209776 455.2066893495879 441.86130196707813 428.5159145845681 415.1362659323852 402.03216741173696 388.8735070919852 375.8033254034093 364.06951641998205 352.3048690470932 342.7985708064939 334.1994146284206 328.11126446807134 323.36200917883906 321.85158340806765 320.3411576372963 321.85158340806765 323.36200917883906 328.11126446807134 334.1994146284206 342.7985708064939 352.3048690470932 364.06951641998205 375.8033254034093 388.8735070919852 402.03216741173696 415.1362659323852 428.5159145845681 441.86130196707813 455.2066893495879 468.55207673209776 481.8841018259439 490.79624246528175 491.4209636271727 491.8619904979994 493.01617917024464 500.5111761579894 500.4547776492527 500.3816528542693 500.2611876427905 496.3233469543589 484.76366010020786 471.75749732833503 458.7464922215454 445.72953530402515 432.69379302804185 419.65751863205605 406.44017404365934 393.13929540174865 379.8213764417027 366.4794486555612 353.1494974328214 340.2637498776179 327.433951108558 315.4237870831155 304.3929798456619 293.37207429232427 284.54018037667265 276.3700768839923 268.50848794620987 264.78561568349545 261.06528451412305 259.2551510053611 257.62757550268054 256.0 256.0 256.0 256.0; 320.0 317.44 314.88 312.32 306.75337754109546 293.5067550821909 280.29701543592626 267.479002885723 254.66360516532086 243.94330773020408 233.39066958953697 223.52800076303424 216.009095165142 208.5019522379521 203.04941823150227 199.37362751581014 196.0029068892851 195.67809370378976 195.4677200030079 196.14904657240712 197.2750955158661 198.98959570612294 201.87544896404046 204.76178298549456 207.72718216972802 210.7082160878625 213.68781445388748 219.27907239784295 231.10140427416457 244.32194660550084 257.5773923394222 270.8214701899742 276.56298150428734 263.7367912164241 250.50668206646597 237.29119381397783 229.61559122417825 229.51984526885036 229.62135414899646 229.7136170963937 229.8032193127928 230.04298322690897 232.3768169441085 234.66541838006162 237.2710844887893 243.6365565079689 249.93022468635115 259.20226633278014 269.3158205412226 281.189136881772 293.4767992554894 306.7383996277447 320.0 333.2616003722553 346.5232007445106 358.810863118228 370.6841794587774 380.79773366721986 390.06977531364885 396.36344349203114 402.7289155112107 405.3345816199384 407.6231830558915 409.957016773091 410.1967806872072 410.2863829036063 410.3786458510035 410.48015473114964 410.3844087758217 402.7088061860222 389.49331793353406 376.2632087835759 363.43701849571266 369.1785298100258 382.4226076605778 395.67805339449916 408.89859572583543 420.7209276021571 426.3121855461125 429.2917839121375 432.27281783027195 435.23821701450544 438.12455103595954 441.01040429387706 442.7249044841339 443.8509534275929 444.5322799969921 444.3219062962103 443.9970931107149 440.62637248418986 436.9505817684977 431.4980477620479 423.990904834858 416.47199923696576 406.60933041046303 396.0566922697959 385.33639483467914 372.520997114277 359.70298456407374 346.4932449178091 333.24662245890454 327.68 325.12 322.56 320.0]
    weights     =   ones(T, size(cps, 2))
    knots       =   T.(knots_vector(deg, size(cps, 2)))
    curve       =   NurbsCurve(cps, knots, weights)
    body        =   DynamicNurbsBody(curve; thk=0, boundary=true)

    return Simulation((6D, 6D), (0,0), D; U, ν, body, T, mem, ϵ)
end
sim         = TwoDimJellyfish(Float32;)
WaterLily.measure_sdf!(sim.flow.σ,sim.body,WaterLily.time(sim))
Plots.contour(sim.flow.σ',levels=[0], label="SD=0",title="SDF contour") 

In my latest simulation, I have checked the velocity field, and the behaviour seems natural with no concerning stuff such as local peak values in cells.

[weymouth]

cps needs to be the right type. Here's a small version with your cps vector showing Float32 and Float64 give identical results.

using Plots, ParametricBodies, StaticArrays
D = 2^7; xs = 0:6D; ys = 0:6D
for T in (Float64,Float32)
    cps = SA{T}[256.0 256.0 256.0 256.0 257.62757550268054 259.2551510053611 261.06528451412305 264.78561568349545 268.50848794620987 276.3700768839923 284.54018037667265 293.37207429232427 304.3929798456619 315.4237870831155 327.433951108558 340.2637498776179 353.1494974328214 366.4794486555612 379.8213764417027 393.13929540174865 406.44017404365934 419.65751863205605 432.69379302804185 445.72953530402515 458.7464922215454 471.75749732833503 484.76366010020786 496.3233469543589 500.2611876427905 500.3816528542693 500.4547776492527 500.5111761579894 493.01617917024464 491.8619904979994 491.4209636271727 490.79624246528175 481.8841018259439 468.55207673209776 455.2066893495879 441.86130196707813 428.5159145845681 415.1362659323852 402.03216741173696 388.8735070919852 375.8033254034093 364.06951641998205 352.3048690470932 342.7985708064939 334.1994146284206 328.11126446807134 323.36200917883906 321.85158340806765 320.3411576372963 321.85158340806765 323.36200917883906 328.11126446807134 334.1994146284206 342.7985708064939 352.3048690470932 364.06951641998205 375.8033254034093 388.8735070919852 402.03216741173696 415.1362659323852 428.5159145845681 441.86130196707813 455.2066893495879 468.55207673209776 481.8841018259439 490.79624246528175 491.4209636271727 491.8619904979994 493.01617917024464 500.5111761579894 500.4547776492527 500.3816528542693 500.2611876427905 496.3233469543589 484.76366010020786 471.75749732833503 458.7464922215454 445.72953530402515 432.69379302804185 419.65751863205605 406.44017404365934 393.13929540174865 379.8213764417027 366.4794486555612 353.1494974328214 340.2637498776179 327.433951108558 315.4237870831155 304.3929798456619 293.37207429232427 284.54018037667265 276.3700768839923 268.50848794620987 264.78561568349545 261.06528451412305 259.2551510053611 257.62757550268054 256.0 256.0 256.0 256.0; 320.0 317.44 314.88 312.32 306.75337754109546 293.5067550821909 280.29701543592626 267.479002885723 254.66360516532086 243.94330773020408 233.39066958953697 223.52800076303424 216.009095165142 208.5019522379521 203.04941823150227 199.37362751581014 196.0029068892851 195.67809370378976 195.4677200030079 196.14904657240712 197.2750955158661 198.98959570612294 201.87544896404046 204.76178298549456 207.72718216972802 210.7082160878625 213.68781445388748 219.27907239784295 231.10140427416457 244.32194660550084 257.5773923394222 270.8214701899742 276.56298150428734 263.7367912164241 250.50668206646597 237.29119381397783 229.61559122417825 229.51984526885036 229.62135414899646 229.7136170963937 229.8032193127928 230.04298322690897 232.3768169441085 234.66541838006162 237.2710844887893 243.6365565079689 249.93022468635115 259.20226633278014 269.3158205412226 281.189136881772 293.4767992554894 306.7383996277447 320.0 333.2616003722553 346.5232007445106 358.810863118228 370.6841794587774 380.79773366721986 390.06977531364885 396.36344349203114 402.7289155112107 405.3345816199384 407.6231830558915 409.957016773091 410.1967806872072 410.2863829036063 410.3786458510035 410.48015473114964 410.3844087758217 402.7088061860222 389.49331793353406 376.2632087835759 363.43701849571266 369.1785298100258 382.4226076605778 395.67805339449916 408.89859572583543 420.7209276021571 426.3121855461125 429.2917839121375 432.27281783027195 435.23821701450544 438.12455103595954 441.01040429387706 442.7249044841339 443.8509534275929 444.5322799969921 444.3219062962103 443.9970931107149 440.62637248418986 436.9505817684977 431.4980477620479 423.990904834858 416.47199923696576 406.60933041046303 396.0566922697959 385.33639483467914 372.520997114277 359.70298456407374 346.4932449178091 333.24662245890454 327.68 325.12 322.56 320.0]
    body = DynamicNurbsBody(BSplineCurve(cps,degree=1))
    vals = [sdf(body, SA{T}[x, y], zero(T)) for y in ys, x in xs]
    sdf_plot = Plots.heatmap(xs, ys, vals; color=:viridis, xlims=(0,6D),ylims=(0,6D), aspect_ratio=1, title="Signed Distance Field for Type $T")
    plot!(body.curve,add_cp=false)
    Plots.contour!(xs, ys, vals, color=:pink, ls=:dash, levels=[0])  # Contour where sdf=0
    display(sdf_plot)
end

As a tip to you, avoid rewriting code (like the uniform knots) and passing arguments you don't need. it makes finding bugs much harder.

[weymouth]

Hey @evanvanderweide. I just put up a PR should fix the NURBS locator for degree>1. Can you try it on your examples?

1. There should be no change for linear splines.
2. The sdf accuracy should be much improved for nonlinear splines.
3. The cost shouldn't be too bad for nonlinear splines. It's still O(length(cps)*length(flow.p)) but the coefficient is certainly bigger than the previous method.

I tried with the example above and these points are correct. Let me know if it's true for your other snap shots.

As I said above, you really should either use fewer control points or use a linear spline. But the locator shouldn't break, in any case.

[evanvanderweide]

Hi @weymouth!

I checked this, and indeed, the SDF plots seem good. I varied it for different numbers of control points between 20 - 50, and all seems good.

However, when I attempt to carry out a simulation, my forces immediately become infinite, and the simulation shuts down, which did not occur before. Then again, I adjusted my simulation setup significantly compared to last time by adding a moving grid, added mass approximations and stuff. I have a simpler test case for a moving circle built out of NURBS curves, so I will run that one later to see if the same happens.

Will come back to you, thanks for the help.

Edit: I ran it for the moving circle and this worked perfectly, so I tried again with the jellyfish case and it seems to be working now, I don't see what went wrong earlier (although I was also working on the sim setup this morning). It does indeed increase the computation time (for my laptop, significantly), so I might stick with the linear splines for my research, as it is also working very well.