From 51a539a291181a78e5d7cf48a1e2efdaa08e3695 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Mon, 4 Nov 2013 18:17:01 -0800 Subject: [PATCH] Generalized to any dimension. Tested. --- SimPEG/mesh/TensorMesh.py | 117 ++++++++++++--------------- SimPEG/tests/test_interpolation.py | 122 ++++++++++++++++++++++++++--- SimPEG/utils/interputils.py | 102 +++++++++++++++++++----- 3 files changed, 244 insertions(+), 97 deletions(-) diff --git a/SimPEG/mesh/TensorMesh.py b/SimPEG/mesh/TensorMesh.py index ea3a0200..87b81303 100644 --- a/SimPEG/mesh/TensorMesh.py +++ b/SimPEG/mesh/TensorMesh.py @@ -315,6 +315,45 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts): # --------------- Methods --------------------- + def getTensor(self, locType): + """ Returns a tensor list. + + :param str locType: What tensor (see below) + :rtype: list + :return: list of the tensors that make up the mesh. + + locType can be:: + + 'Ex' -> x-component of field defined on edges + 'Ey' -> y-component of field defined on edges + 'Ez' -> z-component of field defined on edges + 'Fx' -> x-component of field defined on faces + 'Fy' -> y-component of field defined on faces + 'Fz' -> z-component of field defined on faces + 'N' -> scalar field defined on nodes + 'CC' -> scalar field defined on cell centers + """ + + if locType is 'Fx': + ten = [self.vectorNx , self.vectorCCy, self.vectorCCz] + elif locType is 'Fy': + ten = [self.vectorCCx, self.vectorNy , self.vectorCCz] + elif locType is 'Fz': + ten = [self.vectorCCx, self.vectorCCy, self.vectorNz ] + elif locType is 'Ex': + ten = [self.vectorCCx, self.vectorNy , self.vectorNz ] + elif locType is 'Ey': + ten = [self.vectorNx , self.vectorCCy, self.vectorNz ] + elif locType is 'Ez': + ten = [self.vectorNx , self.vectorNy , self.vectorCCz] + elif locType is 'CC': + ten = [self.vectorCCx, self.vectorCCy, self.vectorCCz] + elif locType is 'N': + ten = [self.vectorNx , self.vectorNy , self.vectorNz ] + + return [t for t in ten if t is not None] + + def isInside(self, pts): """ Determines if a set of points are inside a mesh. @@ -345,9 +384,9 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts): 'Ex' -> x-component of field defined on edges 'Ey' -> y-component of field defined on edges 'Ez' -> z-component of field defined on edges - 'Fx' -> x-component of field defined on edges - 'Fy' -> y-component of field defined on edges - 'Fz' -> z-component of field defined on edges + 'Fx' -> x-component of field defined on faces + 'Fy' -> y-component of field defined on faces + 'Fz' -> z-component of field defined on faces 'N' -> scalar field defined on nodes 'CC' -> scalar field defined on cell centers """ @@ -355,70 +394,16 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts): loc = np.atleast_2d(loc) assert np.all(self.isInside(loc)), "Points outside of mesh" - if self.dim == 3: - - if locType == 'Fx': - Qx = interpmat(self.vectorNx, - self.vectorCCy, - self.vectorCCz, - loc[:,0], loc[:,1], loc[:,2]) - Qy = spzeros(loc.shape[0], self.nF[1]) - Qz = spzeros(loc.shape[0], self.nF[2]) - Q = sp.hstack([Qx, Qy, Qz]) - elif locType == 'Fy': - Qx = spzeros(loc.shape[0], self.nF[0]) - Qy = interpmat(self.vectorCCx, - self.vectorNy, - self.vectorCCz, - loc[:,0], loc[:,1], loc[:,2]) - Qz = spzeros(loc.shape[0], self.nF[2]) - Q = sp.hstack([Qx, Qy, Qz]) - elif locType == 'Fz': - Qx = spzeros(loc.shape[0], self.nF[0]) - Qy = spzeros(loc.shape[0], self.nF[1]) - Qz = interpmat(self.vectorCCx, - self.vectorCCy, - self.vectorNz, - loc[:,0], loc[:,1], loc[:,2]) - Q = sp.hstack([Qx, Qy, Qz]) - elif locType == 'Ex': - Qx = interpmat(self.vectorCCx, - self.vectorNy, - self.vectorNz, - loc[:,0], loc[:,1], loc[:,2]) - Qy = spzeros(loc.shape[0], self.nE[1]) - Qz = spzeros(loc.shape[0], self.nE[2]) - Q = sp.hstack([Qx, Qy, Qz]) - elif locType == 'Ey': - Qx = spzeros(loc.shape[0], self.nE[0]) - Qy = interpmat(self.vectorNx, - self.vectorCCy, - self.vectorNz, - loc[:,0], loc[:,1], loc[:,2]) - Qz = spzeros(loc.shape[0], self.nE[2]) - Q = sp.hstack([Qx, Qy, Qz]) - elif locType == 'Ez': - Qx = spzeros(loc.shape[0], self.nE[0]) - Qy = spzeros(loc.shape[0], self.nE[1]) - Qz = interpmat(self.vectorNx, - self.vectorNy, - self.vectorCCz, - loc[:,0], loc[:,1], loc[:,2]) - Q = sp.hstack([Qx, Qy, Qz]) - elif locType == 'N': - Q = interpmat(self.vectorNx, - self.vectorNy, - self.vectorNz, - loc[:,0], loc[:,1], loc[:,2]) - elif locType == 'CC': - Q = interpmat(self.vectorCCx, - self.vectorCCy, - self.vectorCCz, - loc[:,0], loc[:,1], loc[:,2]) - else: - raise NotImplementedError('getInterpolationMat: locType=='+locType) + ind = 0 if 'x' in locType else 1 if 'y' in locType else 2 if 'z' in locType else -1 + if locType in ['Fx','Fy','Fz','Ex','Ey','Ez'] and self.dim >= ind: + nF_nE = self.nF if 'F' in locType else self.nE + components = [spzeros(loc.shape[0], n) for n in nF_nE] + components[ind] = interpmat(loc, *self.getTensor(locType)) + Q = sp.hstack(components) + elif locType in ['CC', 'N']: + Q = interpmat(loc, *self.getTensor(locType)) else: - raise NotImplementedError('getInterpolationMat: dim=='+str(m.dim)) + raise NotImplementedError('getInterpolationMat: locType=='+locType+' and mesh.dim=='+str(self.dim)) return Q if __name__ == '__main__': diff --git a/SimPEG/tests/test_interpolation.py b/SimPEG/tests/test_interpolation.py index a068c441..6931aa42 100644 --- a/SimPEG/tests/test_interpolation.py +++ b/SimPEG/tests/test_interpolation.py @@ -1,9 +1,11 @@ import numpy as np import unittest from TestUtils import OrderTest +from SimPEG.utils import mkvc MESHTYPES = ['uniformTensorMesh', 'randomTensorMesh'] -TOLERANCES = [0.9, 0.6] +TOLERANCES = [0.9, 0.55] +call1 = lambda fun, xyz: fun(xyz) call2 = lambda fun, xyz: fun(xyz[:, 0], xyz[:, 1]) call3 = lambda fun, xyz: fun(xyz[:, 0], xyz[:, 1], xyz[:, 2]) cart_row2 = lambda g, xfun, yfun: np.c_[call2(xfun, g), call2(yfun, g)] @@ -14,14 +16,114 @@ cartF3 = lambda M, fx, fy, fz: np.vstack((cart_row3(M.gridFx, fx, fy, fz), cart_ cartE3 = lambda M, ex, ey, ez: np.vstack((cart_row3(M.gridEx, ex, ey, ez), cart_row3(M.gridEy, ex, ey, ez), cart_row3(M.gridEz, ex, ey, ez))) -LOCS = np.random.rand(50,3)*0.6+0.2 -class TestInterpolation(OrderTest): +class TestInterpolation1D(OrderTest): + LOCS = np.random.rand(50,1)*0.6+0.2 + name = "Interpolation 1D" + meshTypes = MESHTYPES + tolerance = TOLERANCES + meshDimension = 1 + meshSizes = [8, 16, 32] + + def getError(self): + funX = lambda x: np.cos(2*np.pi*x) + + anal = mkvc(call1(funX, self.LOCS)) + + if 'CC' == self.type: + grid = call1(funX, self.M.gridCC) + elif 'N' == self.type: + grid = call1(funX, self.M.gridN) + + comp = self.M.getInterpolationMat(self.LOCS, self.type)*grid + + err = np.linalg.norm((comp - anal), 2) + return err + + def test_orderCC(self): + self.type = 'CC' + self.name = 'Interpolation 1D: CC' + self.orderTest() + + def test_orderN(self): + self.type = 'N' + self.name = 'Interpolation 1D: N' + self.orderTest() + +class TestInterpolation2d(OrderTest): + name = "Interpolation 2D" + LOCS = np.random.rand(50,2)*0.6+0.2 + meshTypes = MESHTYPES + tolerance = TOLERANCES + meshDimension = 2 + meshSizes = [8, 16, 32, 64] + + def getError(self): + funX = lambda x, y: np.cos(2*np.pi*y) + funY = lambda x, y: np.cos(2*np.pi*x) + + if 'x' in self.type: + anal = call2(funX, self.LOCS) + elif 'y' in self.type: + anal = call2(funY, self.LOCS) + else: + anal = call2(funX, self.LOCS) + + if 'F' in self.type: + Fc = cartF2(self.M, funX, funY) + grid = self.M.projectFaceVector(Fc) + elif 'E' in self.type: + Ec = cartE2(self.M, funX, funY) + grid = self.M.projectEdgeVector(Ec) + elif 'CC' == self.type: + grid = call2(funX, self.M.gridCC) + elif 'N' == self.type: + grid = call2(funX, self.M.gridN) + + comp = self.M.getInterpolationMat(self.LOCS, self.type)*grid + + err = np.linalg.norm((comp - anal), np.inf) + return err + + def test_orderCC(self): + self.type = 'CC' + self.name = 'Interpolation 2D: CC' + self.orderTest() + + def test_orderN(self): + self.type = 'N' + self.name = 'Interpolation 2D: N' + self.orderTest() + + def test_orderFx(self): + self.type = 'Fx' + self.name = 'Interpolation 2D: Fx' + self.orderTest() + + def test_orderFy(self): + self.type = 'Fy' + self.name = 'Interpolation 2D: Fy' + self.orderTest() + + def test_orderEx(self): + self.type = 'Ex' + self.name = 'Interpolation 2D: Ex' + self.orderTest() + + def test_orderEy(self): + self.type = 'Ey' + self.name = 'Interpolation 2D: Ey' + self.orderTest() + + + +class TestInterpolation3D(OrderTest): name = "Interpolation" + LOCS = np.random.rand(50,3)*0.6+0.2 meshTypes = MESHTYPES tolerance = TOLERANCES meshDimension = 3 - meshSizes = [8, 16, 32] + meshSizes = [8, 16, 32, 64] def getError(self): funX = lambda x, y, z: np.cos(2*np.pi*y) @@ -29,13 +131,13 @@ class TestInterpolation(OrderTest): funZ = lambda x, y, z: np.cos(2*np.pi*x) if 'x' in self.type: - anal = call3(funX, LOCS) + anal = call3(funX, self.LOCS) elif 'y' in self.type: - anal = call3(funY, LOCS) + anal = call3(funY, self.LOCS) elif 'z' in self.type: - anal = call3(funZ, LOCS) + anal = call3(funZ, self.LOCS) else: - anal = call3(funX, LOCS) + anal = call3(funX, self.LOCS) if 'F' in self.type: Fc = cartF3(self.M, funX, funY, funZ) @@ -48,7 +150,7 @@ class TestInterpolation(OrderTest): elif 'N' == self.type: grid = call3(funX, self.M.gridN) - comp = self.M.getInterpolationMat(LOCS, self.type)*grid + comp = self.M.getInterpolationMat(self.LOCS, self.type)*grid err = np.linalg.norm((comp - anal), np.inf) return err @@ -94,7 +196,5 @@ class TestInterpolation(OrderTest): self.orderTest() - - if __name__ == '__main__': unittest.main() diff --git a/SimPEG/utils/interputils.py b/SimPEG/utils/interputils.py index 0a544532..0ece29c7 100644 --- a/SimPEG/utils/interputils.py +++ b/SimPEG/utils/interputils.py @@ -3,35 +3,97 @@ import scipy.sparse as sp from sputils import spzeros from matutils import mkvc, sub2ind -def interpmat(x,y,z,xr,yr,zr): +def _interp_point_1D(x, xr_i): + im = np.argmin(abs(x-xr_i)) + if xr_i - x[im] >= 0: # Point on the left + ind_x1 = im + ind_x2 = im+1 + elif xr_i - x[im] < 0: # Point on the right + ind_x1 = im-1 + ind_x2 = im + dx1 = xr_i - x[ind_x1] + dx2 = x[ind_x2] - xr_i + return ind_x1, ind_x2, dx1, dx2 + +def interpmat(locs, x, y=None, z=None): """ Local interpolation computed for each receiver point in turn """ + if y is None and z is None: + return interpmat1D(locs, x) + elif z is None: + return interpmat2D(locs, x, y) + else: + return interpmat3D(locs, x, y, z) + +def interpmat1D(locs, x): + nx = x.size + locs = mkvc(locs) + npts = locs.shape[0] + + Q = sp.lil_matrix((npts, nx)) + + for i in range(npts): + ind_x1, ind_x2, dx1, dx2 = _interp_point_1D(x, locs[i]) + dv = (x[ind_x2] - x[ind_x1]) + Dx = x[ind_x2] - x[ind_x1] + # Get the row in the matrix + inds = [ind_x1, ind_x2] + vals = [(1-dx1/Dx),(1-dx2/Dx)] + Q[i, inds] = vals + return Q.tocsr() + + + +def interpmat2D(locs, x, y): + nx = x.size + ny = y.size + npts = locs.shape[0] + + Q = sp.lil_matrix((npts, nx*ny)) + + + for i in range(npts): + ind_x1, ind_x2, dx1, dx2 = _interp_point_1D(x, locs[i, 0]) + ind_y1, ind_y2, dy1, dy2 = _interp_point_1D(y, locs[i, 1]) + + dv = (x[ind_x2] - x[ind_x1]) * (y[ind_y2] - y[ind_y1]) + + Dx = x[ind_x2] - x[ind_x1] + Dy = y[ind_y2] - y[ind_y1] + + # Get the row in the matrix + + inds = sub2ind((nx,ny),[ + ( ind_x1, ind_y2), + ( ind_x1, ind_y1), + ( ind_x2, ind_y1), + ( ind_x2, ind_y2)]) + + vals = [(1-dx1/Dx)*(1-dy2/Dy), + (1-dx1/Dx)*(1-dy1/Dy), + (1-dx2/Dx)*(1-dy1/Dy), + (1-dx2/Dx)*(1-dy2/Dy)] + + Q[i, mkvc(inds)] = vals + + return Q.tocsr() + + + +def interpmat3D(locs, x, y, z): nx = x.size ny = y.size nz = z.size - npts = xr.shape[0] + npts = locs.shape[0] Q = sp.lil_matrix((npts, nx*ny*nz)) - def inter1D(x, xr_i): - im = np.argmin(abs(x-xr_i)) - if xr_i - x[im] >= 0: # Point on the left - ind_x1 = im - ind_x2 = im+1 - elif xr_i - x[im] < 0: # Point on the right - ind_x1 = im-1 - ind_x2 = im - dx1 = xr_i - x[ind_x1] - dx2 = x[ind_x2] - xr_i - return ind_x1, ind_x2, dx1, dx2 - for i in range(npts): - # in x-direction - ind_x1, ind_x2, dx1, dx2 = inter1D(x, xr[i]) - ind_y1, ind_y2, dy1, dy2 = inter1D(y, yr[i]) - ind_z1, ind_z2, dz1, dz2 = inter1D(z, zr[i]) + ind_x1, ind_x2, dx1, dx2 = _interp_point_1D(x, locs[i, 0]) + ind_y1, ind_y2, dy1, dy2 = _interp_point_1D(y, locs[i, 1]) + ind_z1, ind_z2, dz1, dz2 = _interp_point_1D(z, locs[i, 2]) dv = (x[ind_x2] - x[ind_x1]) * (y[ind_y2] - y[ind_y1]) *(z[ind_z2] - z[ind_z1]) @@ -61,5 +123,5 @@ def interpmat(x,y,z,xr,yr,zr): (1-dx2/Dx)*(1-dy2/Dy)*(1-dz2/Dz)] Q[i, mkvc(inds)] = vals - Q = Q.tocsr() - return Q + + return Q.tocsr()