GOAL

我的目标是使用 python 插值 3D 矢量场。

CODE

原始向量场

import numpy as np
import matplotlib.pyplot as plt
# For interpolation
from scipy.interpolate import RegularGridInterpolator


#%% VECTOR FIELD

xx, yy, zz = np.meshgrid(np.arange(-0.8, 1, 0.2),
                      np.arange(-0.8, 1, 0.2),
                      np.arange(-0.8, 1, 0.8))

uu = np.sin(np.pi * xx) * np.cos(np.pi * yy) * np.cos(np.pi * zz)
vv = -np.cos(np.pi * xx) * np.sin(np.pi * yy) * np.cos(np.pi * zz)
ww = (np.sqrt(2.0 / 3.0) * np.cos(np.pi * xx) * np.cos(np.pi * yy) *
     np.sin(np.pi * zz))

# Ravel -> make 1D lists

x = np.ravel(xx)
y = np.ravel(yy)
z = np.ravel(zz)

u = np.ravel(uu)
v = np.ravel(vv)
w = np.ravel(ww)

插值函数（这里有问题）

#%% INTERPOLATION FUNCTION

def interpolate_field(x,y,z,u,v,w,new_points):
    
    x = np.unique(x)
    y = np.unique(y)
    z = np.unique(z)
    
    u = np.reshape(u, (len(x), len(y), len(z)))
    v = np.reshape(u, (len(x), len(y), len(z)))
    w = np.reshape(u, (len(x), len(y), len(z)))
    
    u_int_f = RegularGridInterpolator((x, y, z), u)
    v_int_f = RegularGridInterpolator((x, y, z), v)
    w_int_f = RegularGridInterpolator((x, y, z), w)
    

    u_int = u_int_f(new_points)
    v_int = v_int_f(new_points)
    w_int = w_int_f(new_points)
    
    return u_int, v_int, w_int

评估新点的插值

#%% EVALUATE INTERPOLATION FUNCTION

new_grid = np.meshgrid(
    np.linspace(np.min(x), np.max(x), 20),
    np.linspace(np.min(y), np.max(y), 20),
    np.linspace(np.min(z), np.max(z), 3)
    , indexing="xy")
            
# create list of new_points
new_points = np.vstack(list(map(np.ravel, new_grid))).T
# get vector field values at new points
uint, vint, wint = interpolate_field(x,y,z,u,v,w,new_points)

new_points = np.array(new_points)
xn = new_points[:,0]
yn = new_points[:,1]
zn = new_points[:,2]

# PLOT
fig = plt.figure(dpi=300)
ax = fig.gca(projection='3d')
ax.quiver(xn, yn, zn, uint, vint, wint, length=0.1)
plt.show()

问题

正如您所看到的，插值函数存在问题，因为生成的矢量图根本没有显示与原始矢量图相同的行为。

这是一种使用的方法RegularGridInterpolator:

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import RegularGridInterpolator


def interp_field(field, coords, shape):
    interpolator = RegularGridInterpolator(
        (gridZ[:, 0, 0], gridY[0, :, 0], gridX[0, 0, :]), field)
    return interpolator(coords).reshape(shape)


gridZ, gridY, gridX = np.mgrid[-0.8:1:3j, -0.8:1:10j, -0.8:1:10j]

u = np.sin(np.pi * gridX) * np.cos(np.pi * gridY) * np.cos(np.pi * gridZ)
v = -np.cos(np.pi * gridX) * np.sin(np.pi * gridY) * np.cos(np.pi * gridZ)
w = (np.sqrt(2.0 / 3.0) * np.cos(np.pi * gridX) * np.cos(np.pi * gridY) *
     np.sin(np.pi * gridZ))

interp_gridZ, interp_gridY, interp_gridX = np.mgrid[
    -0.8:1:3j, -0.8:1:20j, -0.8:1:20j]
interp_coords = np.column_stack(
    (interp_gridZ.flatten(), interp_gridY.flatten(), interp_gridX.flatten()))
interp_u = interp_field(u, interp_coords, interp_gridX.shape)
interp_v = interp_field(v, interp_coords, interp_gridX.shape)
interp_w = interp_field(w, interp_coords, interp_gridX.shape)

fig, (ax, bx) = plt.subplots(ncols=2, subplot_kw=dict(projection='3d'),
                             constrained_layout=True)
ax.quiver(gridX, gridY, gridZ, u, v, w, length=0.3)
bx.quiver(interp_gridX, interp_gridY, interp_gridZ,
          interp_u, interp_v, interp_w, length=0.3)
plt.show()

The resulting plot looks like that: