CODEWITHSUNDEEP
pythonvectorscipy3dinterpolation
henry
henry

Reputation: 975

Correctly interpolate 3D vector field with python scipy

GOAL

My goal is to interpolate a 3D vector field using python.

CODE

Original Vector field

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)

enter image description here

Interpolation Function (something is wrong here)

#%% 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 at new points

#%% 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()

enter image description here

THE ISSUE

As you can see something is wrong with the interpolation function since the resulting vector plot does not show the same behavior as the original at all.

Upvotes: 1

Views: 639

Answers (1)

Patol75
Patol75

Reputation: 4547

Here is a way to do it using 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: enter image description here

Upvotes: 3

Related Questions