用Python插值三维矢量场

Question

任务

我想用python插值一个3D向量场

初始网格->；无关联的u、v、w =>；新网格-->；计算相关的u，v，w

到目前为止我所做的：

创建样本向量场的代码：

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

# REDUCED EXAMPLE
Nx = 8
Ny = 8
Nz = 3
grid = np.meshgrid(np.linspace(-1, 1, Nx),
                   np.linspace(-1, 1, Ny),
                   np.linspace(-1, 1, Nz))

x,y,z = grid

# CONSTRUCT u,v,w 
u = np.sin(np.pi * x) * np.cos(np.pi * y) * np.cos(np.pi * z)
v = -np.cos(np.pi * x) * np.sin(np.pi * y) * np.cos(np.pi * z)
w = (np.sqrt(2.0 / 3.0) * np.cos(np.pi * x) * np.cos(np.pi * y) *
     np.sin(np.pi * z))

fig = plt.figure(dpi=300)
ax = fig.gca(projection='3d')
ax.quiver(x, y, z, u, v, w, length=0.2)
plt.title("Reduced")
plt.show()

。。。我尝试创建一个插值函数：

def interpolate_field(old_points,u,v,w,new_points):
    # old points zip(x,y,z) where x,y,z linspace -> regular grid
    # u,v,w each of shape (NX, NY, NZ)
    # mew points np.array(list(zip(xx,yy,zz))) where xx,yy,zz np.meshgrid
    
    x, y, z = old_points

    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

…但当我尝试应用它时：+

# New grid
NNx = 20
NNy = 20
NNz = 3
grid = np.meshgrid(np.linspace(-1, 1, NNx),
                   np.linspace(-1, 1, NNy),
                   np.linspace(-1, 1, NNz))
    
# Reshape for interpolation function
u_reshape = np.reshape(u, (Nx, Ny, Nz))
v_reshape = np.reshape(v, (Nx, Ny, Nz))
w_reshape = np.reshape(w, (Nx, Ny, Nz))


old_points = (x,y,z)
new_points = np.vstack(list(map(np.ravel, grid))).T

u_int, v_int, w_int = interpolate_field(old_points,
                  u_reshape,
                  v_reshape,
                  w_reshape,
                  new_points)

我得到一个错误：

ValueError: The points in dimension 0 must be strictly ascending

Answer 1

RegularGridInterpolator在应用网格网格之前需要一组点

Parameters points

tuple of ndarray of float, with shapes (m1, ), …, (mn, ) The points defining the regular grid in n dimensions

因此，您可以事先通过网格：

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

def interpolate_field(old_points,u,v,w,new_points):
    # old points zip(x,y,z) where x,y,z linspace -> regular grid
    # u,v,w each of shape (NX, NY, NZ)
    # mew points np.array(list(zip(xx,yy,zz))) where xx,yy,zz np.meshgrid
    
    x, y, z = old_points


    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

# New grid
Nx = 20
Ny = 20
Nz = 3

x = np.linspace(-1, 1, Nx)
y = np.linspace(-1, 1, Ny)
z = np.linspace(-1, 1, Nz)

grid = np.meshgrid(x, y, z)

XX, YY, ZZ = grid

# CONSTRUCT u,v,w 
u = np.sin(np.pi * XX) * np.cos(np.pi * YY) * np.cos(np.pi * ZZ)
v = -np.cos(np.pi * XX) * np.sin(np.pi * YY) * np.cos(np.pi * ZZ)
w = (np.sqrt(2.0 / 3.0) * np.cos(np.pi * XX) * np.cos(np.pi * YY) *
     np.sin(np.pi * ZZ))


# Reshape for interpolation function
u_reshape = np.reshape(u, (Nx, Ny, Nz))
v_reshape = np.reshape(v, (Nx, Ny, Nz))
w_reshape = np.reshape(w, (Nx, Ny, Nz))


old_points = (x,y,z)
new_points = np.vstack(list(map(np.ravel, grid))).T

u_int, v_int, w_int = interpolate_field(old_points,
                  u_reshape,
                  v_reshape,
                  w_reshape,
                  new_points)

或者，您可以尝试从插值_场函数中反转网格网格

编辑：

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

def interpolate_field(old_points,u,v,w,new_points):
    # old points zip(x,y,z) where x,y,z linspace -> regular grid
    # u,v,w each of shape (NX, NY, NZ)
    # mew points np.array(list(zip(xx,yy,zz))) where xx,yy,zz np.meshgrid
    
    x, y, z = old_points
    x = x[0, :, 0]
    y = y[:, 0, 0]
    z = z[0, 0, :]

    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

# New grid
NNx = 20
NNy = 20
NNz = 3

x = np.linspace(-1, 1, NNx)
y = np.linspace(-1, 1, NNy)
z = np.linspace(-1, 1, NNz)

grid = np.meshgrid(np.linspace(-1, 1, NNx),
                   np.linspace(-1, 1, NNy),
                   np.linspace(-1, 1, NNz))

x,y,z = grid

# CONSTRUCT u,v,w 
u = np.sin(np.pi * x) * np.cos(np.pi * y) * np.cos(np.pi * z)
v = -np.cos(np.pi * x) * np.sin(np.pi * y) * np.cos(np.pi * z)
w = (np.sqrt(2.0 / 3.0) * np.cos(np.pi * x) * np.cos(np.pi * y) *
     np.sin(np.pi * z))


# Reshape for interpolation function
u_reshape = np.reshape(u, (NNx, NNy, NNz))
v_reshape = np.reshape(v, (NNx, NNy, NNz))
w_reshape = np.reshape(w, (NNx, NNy, NNz))


old_points = (x,y,z)
new_points = np.vstack(list(map(np.ravel, grid))).T

u_int, v_int, w_int = interpolate_field(old_points,
                  u_reshape,
                  v_reshape,
                  w_reshape,
                    new_points)

u_int = np.reshape(u_int, (NNx,NNy,NNz))
v_int = np.reshape(v_int, (NNx,NNy,NNz))
w_int = np.reshape(w_int, (NNx,NNy,NNz))
fig = plt.figure(dpi=300)
ax = fig.gca(projection='3d')
ax.quiver(x, y, z, u_int, v_int, w_int, length=0.2)
plt.title("Reduced")
plt.show()