# fd_interpolate

## `fd_interpolate(P, gs, h, corner=np.array([0.0, 0.0]))`

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

Returns:

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

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

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``` ``````def fd_interpolate(P,gs,h,corner=np.array([0.0,0.0])): """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)] 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 -------- TO-DO """ dim = P.shape[1] 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