# fd_interpolate

## `fd_interpolate(P, gs, h, corner=None)`

Bi/Trilinear interpolation matrix

Given a regular finite-difference grid described by the number of nodes on each side, the grid spacing, and the location of the bottom-left-front-most corner node, and a list of points, construct a sparse matrix of bilinear interpolation weights so that P = W @ x

Parameters:

Name Type Description Default
`P` `numpy double array `

Matrix of interpolated point coordinates

required
`gs` `numpy int array`

Grid size [nx,ny(,nz)]

required
`h` `numpy double array`

Spacing between grid points [hx,hy(,hz)]

required
`corner`

Location of the bottom-left-front-most corner node

`None`

Returns:

Name Type Description
`W` `scipy sparse.csr_matrix`

Sparse matrix such that if x are the grid nodes, P = W @ x

regular_square_mesh.

#### Notes

The ordering in the output is consistent with the mesh built in regular_square_mesh

Examples:

``````# Grid parameters
gs = np.array([19,15])
h = 1.0/(gs-1)
# Build a grid
x, y = np.meshgrid(np.linspace(0,1,gs[0]),np.linspace(0,1,gs[1]))
V = np.concatenate((np.reshape(x,(-1, 1)),np.reshape(y,(-1, 1))),axis=1)
# Random set of points
P = np.random.rand(10,2)

# Random grid corner
corner = np.random.rand(2)

# Displace by corner
P = P + np.tile(corner,(P.shape[0],1))
V = V + np.tile(corner,(V.shape[0],1))
# Interpolation matrix
from gpytoolbox import fd_interpolate
W = fd_interpolate(P,gs=gs,h=h,corner=corner)
# Now, P = W @ V
``````
Source code in `src/gpytoolbox/fd_interpolate.py`
 ``` 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107``` ``````def fd_interpolate(P,gs,h,corner=None): """Bi/Trilinear interpolation matrix Given a regular finite-difference grid described by the number of nodes on each side, the grid spacing, and the location of the bottom-left-front-most corner node, and a list of points, construct a sparse matrix of bilinear interpolation weights so that P = W @ x Parameters ---------- P : numpy double array Matrix of interpolated point coordinates gs : numpy int array Grid size [nx,ny(,nz)] h : numpy double array Spacing between grid points [hx,hy(,hz)] corner: numpy double array (optional, default None) Location of the bottom-left-front-most corner node Returns ------- W : scipy sparse.csr_matrix Sparse matrix such that if x are the grid nodes, P = W @ x See Also -------- regular_square_mesh. Notes ----- The ordering in the output is consistent with the mesh built in regular_square_mesh Examples -------- ```python # Grid parameters gs = np.array([19,15]) h = 1.0/(gs-1) # Build a grid x, y = np.meshgrid(np.linspace(0,1,gs[0]),np.linspace(0,1,gs[1])) V = np.concatenate((np.reshape(x,(-1, 1)),np.reshape(y,(-1, 1))),axis=1) # Random set of points P = np.random.rand(10,2) # Random grid corner corner = np.random.rand(2) # Displace by corner P = P + np.tile(corner,(P.shape[0],1)) V = V + np.tile(corner,(V.shape[0],1)) # Interpolation matrix from gpytoolbox import fd_interpolate W = fd_interpolate(P,gs=gs,h=h,corner=corner) # Now, P = W @ V ``` """ dim = P.shape[1] if corner is None: corner = np.zeros(dim) indeces = np.floor( (P - np.tile(corner,(P.shape[0],1)))/np.tile(h,(P.shape[0],1)) ).astype(int) I = linspace(0,P.shape[0]-1,P.shape[0],dtype=int) if dim==2: indeces_vectorized = indeces[:,0] + gs[0]*indeces[:,1] I = np.concatenate((I,I,I,I)) J = np.concatenate(( indeces_vectorized,indeces_vectorized+gs[0],indeces_vectorized+1,indeces_vectorized+1+gs[0] )) else: indeces_vectorized = (indeces[:,1] + gs[1]*indeces[:,2])*gs[0] + indeces[:,0] I = np.concatenate((I,I,I,I,I,I,I,I)) J = np.concatenate(( indeces_vectorized,indeces_vectorized+gs[0],indeces_vectorized+1,indeces_vectorized+1+gs[0], indeces_vectorized+(gs[1]*gs[0]),indeces_vectorized+gs[0]+(gs[1]*gs[0]),indeces_vectorized+1+(gs[1]*gs[0]),indeces_vectorized+1+gs[0]+(gs[1]*gs[0]) )) # Position in the bottom left corner vij = np.tile(corner,(P.shape[0],1)) + indeces*np.tile(h,(P.shape[0],1)) # Normalized position inside cell vij = (P - vij)/np.tile(h,(P.shape[0],1)) # Coefficients wrt to each corner of cell if dim==2: coeff_00 = (1-vij[:,1])*(1-vij[:,0]) # bottom left coeff_01 = (1-vij[:,1])*vij[:,0] # bottom right coeff_10 = vij[:,1]*(1-vij[:,0]) # top left coeff_11 = vij[:,1]*vij[:,0] # top right # concatenate (watch that order is consistent with J!!) vals = np.concatenate((coeff_00,coeff_10,coeff_01,coeff_11)) mat_dim = gs[0]*gs[1] else: coeff_000 = (1-vij[:,1])*(1-vij[:,0])*(1-vij[:,2]) # bottom left coeff_010 = (1-vij[:,1])*vij[:,0]*(1-vij[:,2]) # bottom right coeff_100 = vij[:,1]*(1-vij[:,0])*(1-vij[:,2]) # top left coeff_110 = vij[:,1]*vij[:,0]*(1-vij[:,2]) # top right coeff_001 = (1-vij[:,1])*(1-vij[:,0])*vij[:,2] # bottom left coeff_011 = (1-vij[:,1])*vij[:,0]*vij[:,2] # bottom right coeff_101 = vij[:,1]*(1-vij[:,0])*vij[:,2] # top left coeff_111 = vij[:,1]*vij[:,0]*vij[:,2] # top right # concatenate (watch that order is consistent with J!!) vals = np.concatenate((coeff_000,coeff_100,coeff_010,coeff_110,coeff_001,coeff_101,coeff_011,coeff_111)) mat_dim = gs[0]*gs[1]*gs[2] # Build scipy matrix # to-do: maybe add warning if this happens? I = I[(J>=0)*(J=0)*(J=0)*(J