diff --git a/SimPEG/Mesh/CylMesh.py b/SimPEG/Mesh/CylMesh.py index a891d91c..d3ff51be 100644 --- a/SimPEG/Mesh/CylMesh.py +++ b/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__': diff --git a/SimPEG/Mesh/View.py b/SimPEG/Mesh/View.py index 6663b44f..b5ede21f 100644 --- a/SimPEG/Mesh/View.py +++ b/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 = [] diff --git a/SimPEG/Tests/test_cylMesh.py b/SimPEG/Tests/test_cylMesh.py index 5eaef428..8c03380e 100644 --- a/SimPEG/Tests/test_cylMesh.py +++ b/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']