Initial implementation of interpolation matrix generation for TensorMesh (3D ONLY)

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!')

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

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