I have x_1, \dots, x_m points in \mathbb{R}^n defining a polytope as their convex hull, m > n. Typically n = 2 or 3 and m = 4 or 5.
I would like to

make a more or less equispaced (adjusted when necessary) grid \{ z_j \}_j over this polytope,

for an arbitrary point y inside the convex hull, compute positive weights \sum_i w_i = 1 such that for some set of “nearest” points z_1, \dots, z_k of the grid,
y = \sum_i w_i z_i
I am totally fine with whatever approximation/definition of “equispaced” and "nearest’ that make this problem easy. I am solving this problem a few times for a very long computation, on about 50200 gridpoints, so it does not have to be superfast.
I have been looking at packages in JuliaGeometry but this is outside my area of expertise, so I don’t know where to start. Maybe FEM uses methods like this?
(If my domain was a generalized hypercube, I would just cut it up to little hypercubes and use the tensored interpolating weights. But unfortunately it isn’t.)