Evaluate basis at given points

[IvanBioli]

Hi,

Is there a way to evaluate a basis (or any other nutils.function.Array object) in a given set of points? I know that if I have a sample, I can easily evaluate it using the eval method, but can I create a sample using any set of points?

This is a MWE of what I would like to do:

from nutils import mesh
import numpy as np

domain, geom = mesh.rectilinear([np.linspace(0, 1, 10), np.linspace(0, 1, 10)])
basis = domain.basis("spline", degree=2)

# Evaluation with Bezier sample is easy
bezier = domain.sample('bezier', 2)
evaluations = bezier.eval(basis[1])

# But what if want to evaluate basis[1] at the points (0.5, 0.5) and (0.25, 0.25)?

Is there any easy way to obtain this?

Thank you for your help,
Ivan

[joostvanzwieten]

You can create a sample with the desired points using Topology.locate(). Continuing your example:

smpl = domain.locate(geom, [[0.5, 0.5], [0.25, 0.25]], tol=1e-10)
smpl.eval(basis[1])

[IvanBioli]

Thank you for your previous answer, it was exactly what I needed!

I have a follow-up question about the eval method applied to the basis. Should this method return the value of the basis function on the geometry? I'm encountering unexpected results when using a non-uniform knot vector. Specifically, the plot of the basis does not match my expectations—it does not align with the B-Splines on the same knot vector. I have included a MWE below to illustrate this issue. Am I misunderstanding the relationship between the basis and the geometry/topology?

from nutils import mesh
from nutils.expression_v2 import Namespace
import numpy as np
from scipy.interpolate import BSpline
from matplotlib import pyplot as plt

# Define a non-uniform mesh/knot vector
mesh_ = np.linspace(0, 1, 5)
mesh_[2] += 0.10 # <------- Without this the plot is the same
p = 3
knots = np.concatenate((np.zeros(p), mesh_, np.ones(p))) 
nk = knots.size

# Define geometry and topology
topo, geom = mesh.rectilinear([mesh_])
ns = Namespace()
ns.x = geom

# Define function to plot basis (from tutorial)
bezier = topo.sample('bezier', 32)
nanjoin = lambda array, tri: np.insert(array.take(tri.flat, 0).astype(float),
    slice(tri.shape[1], tri.size, tri.shape[1]), np.nan, axis=0)
sampled_x = nanjoin(bezier.eval('x_0' @ ns), bezier.tri)
def plot_line(func, **arguments):
    plt.plot(sampled_x, nanjoin(bezier.eval(func, **arguments), bezier.tri))
    plt.xlabel('x_0')
    plt.xticks(mesh_)

# Plot basis
basis = topo.basis('spline', degree=p)
plot_line(basis)
plt.show()

# Plot B-Splines
for i in range(basis.ndofs):
        c = np.zeros(nk)
        c[i] = 1
        spl = BSpline(knots, c, p)
        plt.plot(sampled_x, spl(sampled_x))

[joostvanzwieten]

In Nutils a Topology has no notion of geometry. If you want to create a B-spline basis on a non-uniform topology-geometry pair, you have to pass the non-uniform knots to Topology.basis:

basis = topo.basis('spline', degree=p, knotvalues=[mesh_])

Note that the knots should be passed without repetitions. Non-default repetitions should be passed as a separate array (multiplicity per vertex, same length as mesh_) via argument knotmultiplicities.

The plot function you copied from the tutorial is rather outdated. I'll try to update the tutorial soon. Instead, you can use nutils.export.triplot. Example with a diff plot (Nutils minus SciPy):

from nutils import export
bezier = topo.sample('bezier', 32)
x, b = bezier.eval([ns.x, basis])
fig = plt.figure()
ax = fig.add_subplot(1,1,1)
for i in range(len(basis)):
    export.triplot(ax, x, values=b[:,i]-BSpline(knots, np.eye(nk)[i], p)(x[:,0]), tri=bezier.tri)
plt.show()