updates to cylindrical to tensor mesh interpolation, tested.

SimPEG/Mesh/CylMesh.py

@@ -1,7 +1,7 @@
 import numpy as np
 import scipy.sparse as sp
 from scipy.constants import pi
-from SimPEG.Utils import mkvc, ndgrid, sdiag, kron3, speye, ddx, av, avExtrap
+from SimPEG.Utils import mkvc, ndgrid, sdiag, kron3, speye, spzeros, ddx, av, avExtrap
 from TensorMesh import BaseTensorMesh
 from InnerProducts import InnerProducts
 from View import CylView
@@ -330,6 +330,62 @@ class CylMesh(BaseTensorMesh, InnerProducts, CylView):
                 raise NotImplementedError('wrapping in the averaging is not yet implemented')
         return self._aveF2CCV
 
+    def getInterpolationMatCartMesh(self, Mrect, locType='CC'):
+        """
+            Takes a cartesian mesh and returns a projection to translate onto the cartesian grid.
+        """
+
+        assert self.isSymmetric, "Currently we have not taken into account other projections for more complicated CylMeshes"
+
+
+        if locType == 'F':
+            # do this three times for each component
+            X = self.getInterpolationMatCartMesh(Mrect, locType='Fx')
+            Y = self.getInterpolationMatCartMesh(Mrect, locType='Fy')
+            Z = self.getInterpolationMatCartMesh(Mrect, locType='Fz')
+            return sp.vstack((X,Y,Z))
+        if locType == 'E':
+            X = self.getInterpolationMatCartMesh(Mrect, locType='Ex')
+            Y = self.getInterpolationMatCartMesh(Mrect, locType='Ey')
+            Z = spzeros(Mrect.nEz, self.nE)
+            return sp.vstack((X,Y,Z))
+
+        grid = getattr(Mrect, 'grid' + locType)
+        # This is unit circle stuff, 0 to 2*pi, starting at x-axis, rotating counter clockwise in an x-y slice
+        theta = - np.arctan2(grid[:,0] - self.cartesianOrigin[0], grid[:,1] - self.cartesianOrigin[1]) + np.pi/2
+        theta[theta < 0] += np.pi*2.0
+        r = ((grid[:,0] - self.cartesianOrigin[0])**2 + (grid[:,1] - self.cartesianOrigin[1])**2)**0.5
+
+        if locType in ['CC', 'N', 'Fz', 'Ez']:
+            G, proj = np.c_[r, theta, grid[:,2]], np.ones(r.size)
+        else:
+            dotMe = {
+                        'Fx': Mrect.normals[:Mrect.nFx,:],
+                        'Fy': Mrect.normals[Mrect.nFx:(Mrect.nFx+Mrect.nFy),:],
+                        'Fz': Mrect.normals[-Mrect.nFz:,:],
+                        'Ex': Mrect.tangents[:Mrect.nEx,:],
+                        'Ey': Mrect.tangents[Mrect.nEx:(Mrect.nEx+Mrect.nEy),:],
+                        'Ez': Mrect.tangents[-Mrect.nEz:,:],
+                    }[locType]
+            if 'F' in locType:
+                normals = np.c_[np.cos(theta), np.sin(theta), np.zeros(theta.size)]
+                proj = ( normals * dotMe ).sum(axis=1)
+            if 'E' in locType:
+                tangents = np.c_[-np.sin(theta), np.cos(theta), np.zeros(theta.size)]
+                proj = ( tangents * dotMe ).sum(axis=1)
+            G = np.c_[r, theta, grid[:,2]]
+
+        interpType = locType
+        if interpType == 'Fy':
+            interpType = 'Fx'
+        elif interpType == 'Ex':
+            interpType = 'Ey'
+
+        Pc2r = self.getInterpolationMat(G, interpType)
+        Proj = sdiag(proj)
+        return Proj * Pc2r
+
+
 
 if __name__ == '__main__':
 

SimPEG/Mesh/View.py

@@ -109,6 +109,10 @@ class TensorView(object):
                 vc = (self.aveN2CC*v).reshape(self.vnC, order='F')
             elif vType in ['Fx', 'Fy', 'Fz', 'Ex', 'Ey', 'Ez']:
                 aveOp = 'ave' + vType[0] + '2CCV'
+                # n = getattr(self,'vn'+vType[0])
+                # if 'x' in vType: v = np.r_[v,np.zeros(n[1]),np.zeros(n[2])]
+                # if 'y' in vType: v = np.r_[np.zeros(n[0]),v,np.zeros(n[2])]
+                # if 'z' in vType: v = np.r_[np.zeros(n[0]),np.zeros(n[1]),v]
                 v = getattr(self,aveOp)*v # average to cell centers
                 ind_xyz = {'x':0,'y':1,'z':2}[vType[1]]
                 vc = self.r(v.reshape((self.nC,-1),order='F'), 'CC','CC','M')[ind_xyz]
@@ -190,9 +194,16 @@ class TensorView(object):
             M.plotSlice(M.cellGrad*b, 'F', view='vec', grid=True, showIt=True, pcolorOpts={'alpha':0.8})
 
         """
+        if type(vType) in [list, tuple]:
+            assert ax is None, "cannot specify an axis to plot on with this function."
+            fig, axs = plt.subplots(1,len(vType))
+            out = []
+            for vTypeI, ax in zip(vType, axs):
+                out += [self.plotSlice(v,vType=vTypeI, normal=normal, ind=ind, grid=grid, view=view, ax=ax, clim=clim, showIt=False, pcolorOpts=pcolorOpts, streamOpts=streamOpts, gridOpts=gridOpts)]
+            return out
         viewOpts = ['real','imag','abs','vec']
         normalOpts = ['X', 'Y', 'Z']
-        vTypeOpts = ['CC', 'CCv','F','E']
+        vTypeOpts = ['CC', 'CCv','F','E','Fx','Fy','Fz','E','Ex','Ey','Ez']
 
         # Some user error checking
         assert vType in vTypeOpts, "vType must be in ['%s']" % "','".join(vTypeOpts)
@@ -219,9 +230,15 @@ class TensorView(object):
             elif vType == 'CCv':
                 assert view == 'vec', 'Other types for CCv not supported'
             else:
-                # Now just deal with 'F' and 'E'
+                # Now just deal with 'F' and 'E' (x,y,z, maybe...)
                 aveOp = 'ave' + vType + ('2CCV' if view == 'vec' else '2CC')
-                v = getattr(self,aveOp)*v # average to cell centers (might be a vector)
+                Av = getattr(self,aveOp)
+                if v.size == Av.shape[1]:
+                    v = Av * v
+                else:
+                    v = self.r(v,vType[0],vType) # get specific component
+                    v = Av * v
+                # we should now be averaged to cell centers (might be a vector)
             v = self.r(v.reshape((self.nC,-1),order='F'),'CC','CC','M')
             if view == 'vec':
                 outSlice = []

SimPEG/Tests/test_cylMesh.py

@@ -133,6 +133,80 @@ class TestCyl2DMesh(unittest.TestCase):
     def test_lightOperators(self):
         self.assertTrue(self.mesh.nodalGrad is None)
 
+    def test_getInterpMatCartMesh_Cells(self):
+
+        Mr = Mesh.TensorMesh([100,100,2], x0='CC0')
+        Mc = Mesh.CylMesh([np.ones(10)/5,1,10],x0='0C0',cartesianOrigin=[-0.2,-0.2,0])
+
+        mc = np.arange(Mc.nC)
+        xr = np.linspace(0,0.4,50)
+        xc = np.linspace(0,0.4,50) + 0.2
+        Pr = Mr.getInterpolationMat(np.c_[xr,np.ones(50)*-0.2,np.ones(50)*0.5],'CC')
+        Pc = Mc.getInterpolationMat(np.c_[xc,np.zeros(50),np.ones(50)*0.5],'CC')
+        Pc2r = Mc.getInterpolationMatCartMesh(Mr, 'CC')
+
+        assert np.abs(Pr*(Pc2r*mc) - Pc*mc).max() < 1e-3
+
+    def test_getInterpMatCartMesh_Faces(self):
+
+        Mr = Mesh.TensorMesh([100,100,2], x0='CC0')
+        Mc = Mesh.CylMesh([np.ones(10)/5,1,10],x0='0C0',cartesianOrigin=[-0.2,-0.2,0])
+
+        Pf = Mc.getInterpolationMatCartMesh(Mr, 'F')
+        mf = np.ones(Mc.nF)
+
+        frect = Pf * mf
+
+        fxcc = Mr.aveFx2CC*Mr.r(frect, 'F', 'Fx')
+        fycc = Mr.aveFy2CC*Mr.r(frect, 'F', 'Fy')
+        fzcc = Mr.r(frect, 'F', 'Fz')
+
+        indX = Utils.closestPoints(Mr, [0.45, -0.2, 0.5])
+        indY = Utils.closestPoints(Mr, [-0.2, 0.45, 0.5])
+
+        TOL = 1e-2
+        assert np.abs(float(fxcc[indX]) - 1) < TOL
+        assert np.abs(float(fxcc[indY]) - 0) < TOL
+        assert np.abs(float(fycc[indX]) - 0) < TOL
+        assert np.abs(float(fycc[indY]) - 1) < TOL
+        assert np.abs((fzcc - 1).sum()) < TOL
+
+        mag = (fxcc**2 + fycc**2)**0.5
+        dist = ((Mr.gridCC[:,0] + 0.2)**2  + (Mr.gridCC[:,1] + 0.2)**2)**0.5
+
+        assert np.abs(mag[dist > 0.1].max() - 1) < TOL
+        assert np.abs(mag[dist > 0.1].min() - 1) < TOL
+
+
+    def test_getInterpMatCartMesh_Edges(self):
+
+        Mr = Mesh.TensorMesh([100,100,2], x0='CC0')
+        Mc = Mesh.CylMesh([np.ones(10)/5,1,10],x0='0C0',cartesianOrigin=[-0.2,-0.2,0])
+
+        Pe = Mc.getInterpolationMatCartMesh(Mr, 'E')
+        me = np.ones(Mc.nE)
+
+        erect = Pe * me
+
+        excc = Mr.aveEx2CC*Mr.r(erect, 'E', 'Ex')
+        eycc = Mr.aveEy2CC*Mr.r(erect, 'E', 'Ey')
+        ezcc = Mr.r(erect, 'E', 'Ez')
+
+        indX = Utils.closestPoints(Mr, [0.45, -0.2, 0.5])
+        indY = Utils.closestPoints(Mr, [-0.2, 0.45, 0.5])
+
+        TOL = 1e-2
+        assert np.abs(float(excc[indX]) - 0) < TOL
+        assert np.abs(float(excc[indY]) + 1) < TOL
+        assert np.abs(float(eycc[indX]) - 1) < TOL
+        assert np.abs(float(eycc[indY]) - 0) < TOL
+        assert np.abs(ezcc.sum()) < TOL
+
+        mag = (excc**2 + eycc**2)**0.5
+        dist = ((Mr.gridCC[:,0] + 0.2)**2  + (Mr.gridCC[:,1] + 0.2)**2)**0.5
+
+        assert np.abs(mag[dist > 0.1].max() - 1) < TOL
+        assert np.abs(mag[dist > 0.1].min() - 1) < TOL
 
 
 MESHTYPES = ['uniformCylMesh']