numpy array export

Is there a built-in function for interpolating a GridFunction/CoefficientFunction onto an n-dimensional grid for export to numpy, i.e. the inverse of VoxelCoefficient? There are a lot of nice built-in functions for evaluation of data, but there are limits, and I think at some point I will need to use numpy for analysis. I've searched tutorials and the forum and only found one user-written function to do this that I'm not completely sure about.

I ended up adapting a function posted by s.mokbel here for 3D with a scalar u, and a touch more robustness (np.nan for outside mesh, and uses maximum and minimum points in the mesh for the interpolation).

def create_np_data(gfu: GridFunction, mesh: Mesh, l: int, m: int, n: int):
    """
    Function to create numpy array from Gridfunction data.
    Args:
    gfu: Scalar gridfunction
    mesh: The mesh used in this simulation.
    l: Desired number of discrete cells (elements) in the uniform grid (numpy array)
    in the x direction.
    n: Desired number of discrete cells (elements) in the uniform grid (numpy array)
    in the y direction.
    m: Desired number of discrete cells (elements) in the uniform grid (numpy array)
    in the z direction.
    Returns:
    Numpy array of field data U
    """
    # Obtain mesh dimensions from vertices.
    vertices = np.array([vertex.point for vertex in mesh.vertices])
    max_vertex = vertices.max(axis=0)
    min_vertex = vertices.min(axis=0)
    x0, y0, z0 = min_vertex
    xn, yn, zn = max_vertex
    # Initialize numpy arrays for field data.
    output_U = np.zeros((l,m,n),dtype=float)
    # Create uniform grid points based on mesh dimensions.
    x_interp = np.linspace(x0,xn,l)
    y_interp = np.linspace(y0,yn,m)
    z_interp = np.linspace(z0,zn,n)
                           
    x_idx = 0
    for p1 in x_interp:
        y_idx = 0
        for p2 in y_interp:
            z_idx = 0
            for p3 in z_interp:
                try:
                    val_U = gfu(mesh(p1,p2,p3))
                    output_U[x_idx, y_idx, z_idx] = val_U
                except Exception:
                    # Catch any points that might not be in the mesh,
                    # and set their values to np.nan
                    output_U[x_idx, y_idx, z_idx] = np.nan
                z_idx += 1
            y_idx += 1
        x_idx += 1
    return output_U