diff --git a/src/ttcrpy/rgrid.pyx b/src/ttcrpy/rgrid.pyx index 1197715..12f036e 100644 --- a/src/ttcrpy/rgrid.pyx +++ b/src/ttcrpy/rgrid.pyx @@ -577,7 +577,7 @@ cdef class Grid3d: raise ValueError('velocity must be 1D or 3D ndarray') self.grid.setSlowness(slown) - def compute_D(self, coord): + def compute_D(self, np.ndarray[np.double_t, ndim=2] coord): """ compute_D(coord) @@ -593,6 +593,11 @@ cdef class Grid3d: ------- D : scipy csr_matrix, shape (npts, nparams) Matrix of interpolation weights + + Note + ---- + In the current implementation, no check is made to see if the coordinates + are on a node, edge, face, or corner. """ if self.is_outside(coord): raise ValueError('Velocity data point outside grid') @@ -608,8 +613,8 @@ cdef class Grid3d: vec = np.ones(ivec.shape) for n in np.arange(coord.shape[0]): i = int((coord[n,0]-self.x[0])/self._dx ) - j = int((coord[n,1]-self.y[0])/self._dx ) - k = int((coord[n,2]-self.z[0])/self._dx ) + j = int((coord[n,1]-self.y[0])/self._dy ) + k = int((coord[n,2]-self.z[0])/self._dz ) jvec[n] = self.indc(i,j,k) @@ -2565,6 +2570,71 @@ cdef class Grid2d: else: raise ValueError('v must be 1D or 2D ndarray') self.grid.setR4(data) + + def compute_D(self, np.ndarray[np.double_t, ndim=2] coord): + """ + compute_D(coord) + + Return matrix of interpolation weights for velocity data points + constraint + + Parameters + ---------- + coord : np.ndarray, shape (npts, 2) + coordinates of data points + + Returns + ------- + D : scipy csr_matrix, shape (npts, nparams) + Matrix of interpolation weights + + Note + ---- + In the current implementation, no check is made to see if the coordinates + are on a node, edge, or corner. + """ + if self.is_outside(coord): + raise ValueError('Velocity data point outside grid') + + # D is npts x nparams + + cdef Py_ssize_t n, i, k, i1, i2, k1, k2, ii + + if self.cell_slowness: + # for each point in coord, we have 1 values in D + ivec = np.arange(coord.shape[0], dtype=np.int64) + jvec = np.zeros(ivec.shape, dtype=np.int64) + vec = np.ones(ivec.shape) + for n in np.arange(coord.shape[0]): + i = int((coord[n,0]-self.x[0])/self._dx ) + k = int((coord[n,1]-self.z[0])/self._dz ) + + jvec[n] = i * (self._z.size()-1) + k + + return sp.csr_matrix((vec, (ivec,jvec)), + shape=(coord.shape[0], self.get_number_of_cells())) + else: + # for each point in coord, we have 4 values in D + ivec = np.kron(np.arange(coord.shape[0], dtype=np.int64),np.ones(4, dtype=np.int64)) + jvec = np.zeros(ivec.shape, dtype=np.int64) + vec = np.zeros(ivec.shape) + for n in np.arange(coord.shape[0]): + i1 = int(1.e-6 + (coord[n,0]-self._x[0])/self._dx ) + i2 = i1 + 1 + k1 = int(1.e-6 + (coord[n,1]-self._z[0])/self._dz ) + k2 = k1 + 1 + + ii = 0 + for i in (i1, i2): + for k in (k1, k2): + jvec[n*4+ii] = i * self._z.size() + k + vec[n*4+ii] = ((1. - np.abs(coord[n,0]-self._x[i])/self._dx) * + (1. - np.abs(coord[n,1]-self._z[k])/self._dz)) + ii += 1 + + return sp.csr_matrix((vec, (ivec,jvec)), + shape=(coord.shape[0], self.get_number_of_nodes())) + def compute_K(self, order=1): """ diff --git a/tests/test_rgrid2d.py b/tests/test_rgrid2d.py index b84f78a..00a6fe7 100644 --- a/tests/test_rgrid2d.py +++ b/tests/test_rgrid2d.py @@ -216,7 +216,83 @@ def test_Grid2Dsp(self): tt_ref = get_tt('./files/Grid2Drcsp_tt_grid_weakly.vtr') self.assertLess(np.sum(np.abs(tt-tt_ref))/tt.size, 0.01, 'SPM accuracy failed (weakly anelliptical)') +# +class TestComputeD(unittest.TestCase): + + def test_cell_slowness(self): + reader = vtk.vtkXMLRectilinearGridReader() + + reader.SetFileName('./files/layers_fine2d.vtr') + reader.Update() + + data = reader.GetOutput() + x = vtk_to_numpy(data.GetXCoordinates()) + z = vtk_to_numpy(data.GetZCoordinates()) + + slowness = vtk_to_numpy(data.GetCellData().GetArray('Slowness')) + dim = (x.size - 1, z.size - 1) + slowness = slowness.reshape(dim, order='F').flatten() + + g = rg.Grid2d(x, z, method='FSM') + + rng = np.random.default_rng() + coord = np.c_[rng.uniform(0.001, 20.0, 100), rng.uniform(0.001, 20.0, 100)] + D = g.compute_D(coord) + s1 = D @ slowness + + pts = vtk.vtkPoints() + for n in range(coord.shape[0]): + pts.InsertNextPoint(coord[n, 0], 0.0, coord[n, 1]) + + ppts = vtk.vtkPolyData() + ppts.SetPoints(pts) + + pf = vtk.vtkProbeFilter() + pf.SetInputData(ppts) + pf.SetSourceData(data) + pf.Update() + + s2 = vtk_to_numpy(pf.GetOutput().GetPointData().GetArray('Slowness')) + + self.assertAlmostEqual(np.sum(np.abs(s1 - s2)), 0.0) + + def test_node_slowness(self): + reader = vtk.vtkXMLRectilinearGridReader() + + reader.SetFileName('./files/gradient_fine2d.vtr') + reader.Update() + + data = reader.GetOutput() + x = vtk_to_numpy(data.GetXCoordinates()) + z = vtk_to_numpy(data.GetZCoordinates()) + + slowness = vtk_to_numpy(data.GetPointData().GetArray('Slowness')) + dim = (x.size, z.size) + slowness = slowness.reshape(dim, order='F').flatten() + + g = rg.Grid2d(x, z, method='FSM', cell_slowness=0) + + rng = np.random.default_rng() + coord = np.c_[rng.uniform(0.001, 20.0, 100), rng.uniform(0.001, 20.0, 100)] + D = g.compute_D(coord) + s1 = D @ slowness + + pts = vtk.vtkPoints() + for n in range(coord.shape[0]): + pts.InsertNextPoint(coord[n, 0], 0.0, coord[n, 1]) + + ppts = vtk.vtkPolyData() + ppts.SetPoints(pts) + + interpolator = vtk.vtkPointInterpolator() + interpolator.SetInputData(ppts) + interpolator.SetSourceData(data) + interpolator.Update() + + s2 = vtk_to_numpy(interpolator.GetOutput().GetPointData().GetArray('Slowness')) + + self.assertLess(np.sum(np.abs(s1-s2))/coord.shape[0], 0.01, "compute_D accuracy failed slowness at nodes") class Data_kernel(unittest.TestCase):