diff --git a/SimPEG/mesh/TensorMesh.py b/SimPEG/mesh/TensorMesh.py index a3329c14..1dbd62a3 100644 --- a/SimPEG/mesh/TensorMesh.py +++ b/SimPEG/mesh/TensorMesh.py @@ -1,9 +1,10 @@ import numpy as np +import scipy.sparse as sp from BaseMesh import BaseMesh from TensorView import TensorView from DiffOperators import DiffOperators from InnerProducts import InnerProducts -from SimPEG.utils import ndgrid, mkvc +from SimPEG.utils import ndgrid, mkvc, spzeros, interpmat class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts): @@ -312,6 +313,113 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts): _edge = None edge = property(**edge()) + # --------------- Methods --------------------- + + def isInside(self, pts): + """ + Determines if a set of points are inside a mesh. + + :param numpy.ndarray pts: Location of points to test + :rtype numpy.ndarray + :return inside, numpy array of booleans + """ + + pts = np.atleast_2d(pts) + inside = (pts[:,0] >= self.vectorNx.min()) & (pts[:,0] <= self.vectorNx.max()) + if self.dim > 1: + inside = inside & ((pts[:,1] >= self.vectorNy.min()) & (pts[:,1] <= self.vectorNy.max())) + if self.dim > 2: + inside = inside & ((pts[:,2] >= self.vectorNz.min()) & (pts[:,2] <= self.vectorNz.max())) + return inside + + def getInterpolationMat(self, loc, locType): + """ Produces interpolation matrix + + :param numpy.ndarray loc: Location of points to interpolate to + :param str locType: What to interpolate (see below) + :rtype: scipy.sparse.csr.csr_matrix + :return: M, the interpolation matrix + + 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 edges + 'fy' -> y-component of field defined on edges + 'fz' -> z-component of field defined on edges + 'n' -> scalar field defined on nodes + 'cc' -> scalar field defined on cell centres + """ + + 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.nF[1]) + Qz = spzeros(loc.shape[0], self.nF[2]) + Q = sp.hstack([Qx, Qy, Qz]) + elif locType == 'ey': + Qx = spzeros(loc.shape[0], self.nF[0]) + Qy = interpmat(self.vectorNx, + self.vectorCCy, + self.vectorNz, + loc[:,0], loc[:,1], loc[:,2]) + Qz = spzeros(loc.shape[0], self.nF[2]) + Q = sp.hstack([Qx, Qy, Qz]) + elif locType == 'ez': + Qx = spzeros(loc.shape[0], self.nF[0]) + Qy = spzeros(loc.shape[0], self.nF[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) + else: + raise NotImplementedError('getInterpolationMat: dim=='+str(m.dim)) + return Q if __name__ == '__main__': print('Welcome to tensor mesh!') diff --git a/SimPEG/utils/__init__.py b/SimPEG/utils/__init__.py index 7c976ce3..9457df03 100644 --- a/SimPEG/utils/__init__.py +++ b/SimPEG/utils/__init__.py @@ -1,7 +1,9 @@ import matutils import sputils import lomutils +import interputils import ModelBuilder from matutils import getSubArray, mkvc, ndgrid, ind2sub, sub2ind from sputils import spzeros, kron3, speye, sdiag from lomutils import volTetra, faceInfo, inv2X2BlockDiagonal, inv3X3BlockDiagonal, indexCube, exampleLomGird +from interputils import interpmat \ No newline at end of file diff --git a/SimPEG/utils/interputils.py b/SimPEG/utils/interputils.py new file mode 100644 index 00000000..0d528e93 --- /dev/null +++ b/SimPEG/utils/interputils.py @@ -0,0 +1,77 @@ +import numpy as np +import scipy.sparse as sp +from sputils import spzeros +from matutils import mkvc, sub2ind + +def interpmat(x,y,z,xr,yr,zr): + + """ Local interpolation computed for each receiver point in turn """ + + nx = max(x.shape) + ny = max(y.shape) + nz = max(z.shape) + npts = max(xr.shape) + + Q = sp.lil_matrix((npts, nx*ny*nz)) + + for i in range(npts): + # in x-direction + 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] + # in y-direction + im = np.argmin(abs(y-yr[i])) + if yr[i] - y[im] >= 0: # Point on the left + ind_y1 = im + ind_y2 = im+1 + elif yr[i] - y[im] < 0: # Point on the right + ind_y1 = im-1 + ind_y2 = im + dy1 = yr[i] - y[ind_y1] + dy2 = y[ind_y2] - yr[i] + # in z-direction + im = np.argmin(abs(z-zr[i])) + if zr[i] - z[im] >= 0: # Point on the left + ind_z1 = im + ind_z2 = im+1 + elif zr[i] - z[im] < 0: # Point on the right + ind_z1 = im-1 + ind_z2 = im + dz1 = zr[i] - z[ind_z1] + dz2 = z[ind_z2] - zr[i] + dv = (x[ind_x2] - x[ind_x1]) * (y[ind_y2] - y[ind_y1]) *(z[ind_z2] - z[ind_z1]) + + Dx = x[ind_x2] - x[ind_x1] + Dy = y[ind_y2] - y[ind_y1] + Dz = z[ind_z2] - z[ind_z1] + + # Get the row in the matrix + + inds = sub2ind((nx,ny,nz),[ + ( ind_x1, ind_y2, ind_z1), + ( ind_x1, ind_y1, ind_z1), + ( ind_x2, ind_y1, ind_z1), + ( ind_x2, ind_y2, ind_z1), + ( ind_x1, ind_y1, ind_z2), + ( ind_x1, ind_y2, ind_z2), + ( ind_x2, ind_y1, ind_z2), + ( ind_x2, ind_y2, ind_z2)]) + + vals = [(1-dx1/Dx)*(1-dy2/Dy)*(1-dz1/Dz), + (1-dx1/Dx)*(1-dy1/Dy)*(1-dz1/Dz), + (1-dx2/Dx)*(1-dy1/Dy)*(1-dz1/Dz), + (1-dx2/Dx)*(1-dy2/Dy)*(1-dz1/Dz), + (1-dx1/Dx)*(1-dy1/Dy)*(1-dz2/Dz), + (1-dx1/Dx)*(1-dy2/Dy)*(1-dz2/Dz), + (1-dx2/Dx)*(1-dy1/Dy)*(1-dz2/Dz), + (1-dx2/Dx)*(1-dy2/Dy)*(1-dz2/Dz)] + + Q[i, mkvc(inds)] = vals + Q = Q.tocsr() + return Q \ No newline at end of file