From 16c269e01350b856e3e64e36cd200b7530167851 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Thu, 19 Mar 2015 12:48:50 -0700 Subject: [PATCH 1/3] updates to cyl mesh and minor bug fixes. --- SimPEG/Mesh/CylMesh.py | 5 ++++- SimPEG/Mesh/TensorMesh.py | 20 +++++++++++--------- 2 files changed, 15 insertions(+), 10 deletions(-) diff --git a/SimPEG/Mesh/CylMesh.py b/SimPEG/Mesh/CylMesh.py index 8d2b1ae8..a891d91c 100644 --- a/SimPEG/Mesh/CylMesh.py +++ b/SimPEG/Mesh/CylMesh.py @@ -27,11 +27,14 @@ class CylMesh(BaseTensorMesh, InnerProducts, CylView): _unitDimensions = [1, 2*np.pi, 1] - def __init__(self, h, x0=None): + def __init__(self, h, x0=None, cartesianOrigin=None): BaseTensorMesh.__init__(self, h, x0) assert self.hy.sum() == 2*np.pi, "The 2nd dimension must sum to 2*pi" if self.dim == 2: print 'Warning, a disk mesh has not been tested thoroughly.' + cartesianOrigin = np.zeros(self.dim) if cartesianOrigin is None else cartesianOrigin + assert len(cartesianOrigin) == self.dim, "cartesianOrigin must be the same length as the dimension of the mesh." + self.cartesianOrigin = np.array(cartesianOrigin, dtype=float) @property diff --git a/SimPEG/Mesh/TensorMesh.py b/SimPEG/Mesh/TensorMesh.py index cd46fb85..37bfea86 100644 --- a/SimPEG/Mesh/TensorMesh.py +++ b/SimPEG/Mesh/TensorMesh.py @@ -168,21 +168,21 @@ class BaseTensorMesh(BaseRectangularMesh): """ - if key is 'Fx': + if key == 'Fx': ten = [self.vectorNx , self.vectorCCy, self.vectorCCz] - elif key is 'Fy': + elif key == 'Fy': ten = [self.vectorCCx, self.vectorNy , self.vectorCCz] - elif key is 'Fz': + elif key == 'Fz': ten = [self.vectorCCx, self.vectorCCy, self.vectorNz ] - elif key is 'Ex': + elif key == 'Ex': ten = [self.vectorCCx, self.vectorNy , self.vectorNz ] - elif key is 'Ey': + elif key == 'Ey': ten = [self.vectorNx , self.vectorCCy, self.vectorNz ] - elif key is 'Ez': + elif key == 'Ez': ten = [self.vectorNx , self.vectorNy , self.vectorCCz] - elif key is 'CC': + elif key == 'CC': ten = [self.vectorCCx, self.vectorCCy, self.vectorCCz] - elif key is 'N': + elif key == 'N': ten = [self.vectorNx , self.vectorNy , self.vectorNz ] return [t for t in ten if t is not None] @@ -204,10 +204,12 @@ class BaseTensorMesh(BaseRectangularMesh): if locType == 'N' and self._meshType == 'CYL': #NOTE: for a CYL mesh we add a node to check if we are inside in the radial direction! tensors[0] = np.r_[0.,tensors[0]] + tensors[1] = np.r_[tensors[1], 2.0*np.pi] inside = np.ones(pts.shape[0],dtype=bool) for i, tensor in enumerate(tensors): - inside = inside & (pts[:,i] >= tensor.min()) & (pts[:,i] <= tensor.max()) + TOL = np.diff(tensor).min() * 1.0e-10 + inside = inside & (pts[:,i] >= tensor.min()-TOL) & (pts[:,i] <= tensor.max()+TOL) return inside def getInterpolationMat(self, loc, locType, zerosOutside=False): From 354a137461dab5133ac6963c4fdf05041880cb5d Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Thu, 19 Mar 2015 13:18:45 -0700 Subject: [PATCH 2/3] add additional operators to the tensor mesh. --- SimPEG/Mesh/DiffOperators.py | 131 ++++++++++++++++++++++++----------- 1 file changed, 89 insertions(+), 42 deletions(-) diff --git a/SimPEG/Mesh/DiffOperators.py b/SimPEG/Mesh/DiffOperators.py index e46bed20..00bf0259 100644 --- a/SimPEG/Mesh/DiffOperators.py +++ b/SimPEG/Mesh/DiffOperators.py @@ -571,35 +571,58 @@ class DiffOperators(object): @property def aveF2CC(self): "Construct the averaging operator on cell faces to cell centers." - if getattr(self, '_aveF2CC', None) is None: - n = self.vnC - if(self.dim == 1): - self._aveF2CC = av(n[0]) - elif(self.dim == 2): - self._aveF2CC = (0.5)*sp.hstack((sp.kron(speye(n[1]), av(n[0])), - sp.kron(av(n[1]), speye(n[0]))), format="csr") - elif(self.dim == 3): - self._aveF2CC = (1./3.)*sp.hstack((kron3(speye(n[2]), speye(n[1]), av(n[0])), - kron3(speye(n[2]), av(n[1]), speye(n[0])), - kron3(av(n[2]), speye(n[1]), speye(n[0]))), format="csr") - return self._aveF2CC - + if(self.dim == 1): + return self.aveFx2CC + elif(self.dim == 2): + return (0.5)*sp.hstack((self.aveFx2CC, self.aveFy2CC), format="csr") + elif(self.dim == 3): + return (1./3.)*sp.hstack((self.aveFx2CC, self.aveFy2CC, self.aveFz2CC), format="csr") @property def aveF2CCV(self): "Construct the averaging operator on cell faces to cell centers." - if getattr(self, '_aveF2CCV', None) is None: + if(self.dim == 1): + return self.aveFx2CC + elif(self.dim == 2): + return sp.block_diag((self.aveFx2CC, self.aveFy2CC), format="csr") + elif(self.dim == 3): + return sp.block_diag((self.aveFx2CC, self.aveFy2CC, self.aveFz2CC), format="csr") + + @property + def aveFx2CC(self): + "Construct the averaging operator on cell faces in the x direction to cell centers." + if getattr(self, '_aveFx2CC', None) is None: n = self.vnC if(self.dim == 1): - self._aveF2CCV = av(n[0]) + self._aveFx2CC = av(n[0]) elif(self.dim == 2): - self._aveF2CCV = sp.block_diag((sp.kron(speye(n[1]), av(n[0])), - sp.kron(av(n[1]), speye(n[0]))), format="csr") + self._aveFx2CC = sp.kron(speye(n[1]), av(n[0])) elif(self.dim == 3): - self._aveF2CCV = sp.block_diag((kron3(speye(n[2]), speye(n[1]), av(n[0])), - kron3(speye(n[2]), av(n[1]), speye(n[0])), - kron3(av(n[2]), speye(n[1]), speye(n[0]))), format="csr") - return self._aveF2CCV + self._aveFx2CC = kron3(speye(n[2]), speye(n[1]), av(n[0])) + return self._aveFx2CC + + @property + def aveFy2CC(self): + "Construct the averaging operator on cell faces in the y direction to cell centers." + if self.dim < 2: return None + if getattr(self, '_aveFy2CC', None) is None: + n = self.vnC + if(self.dim == 2): + self._aveFy2CC = sp.kron(av(n[1]), speye(n[0])) + elif(self.dim == 3): + self._aveFy2CC = kron3(speye(n[2]), av(n[1]), speye(n[0])) + return self._aveFy2CC + + @property + def aveFz2CC(self): + "Construct the averaging operator on cell faces in the z direction to cell centers." + if self.dim < 3: return None + if getattr(self, '_aveFz2CC', None) is None: + n = self.vnC + if(self.dim == 3): + self._aveFz2CC = kron3(av(n[2]), speye(n[1]), speye(n[0])) + return self._aveFz2CC + @property def aveCC2F(self): @@ -620,36 +643,60 @@ class DiffOperators(object): @property def aveE2CC(self): "Construct the averaging operator on cell edges to cell centers." - if getattr(self, '_aveE2CC', None) is None: - # The number of cell centers in each direction - n = self.vnC - if(self.dim == 1): - self._aveE2CC = speye(n[0]) - elif(self.dim == 2): - self._aveE2CC = 0.5*sp.hstack((sp.kron(av(n[1]), speye(n[0])), - sp.kron(speye(n[1]), av(n[0]))), format="csr") - elif(self.dim == 3): - self._aveE2CC = (1./3)*sp.hstack((kron3(av(n[2]), av(n[1]), speye(n[0])), - kron3(av(n[2]), speye(n[1]), av(n[0])), - kron3(speye(n[2]), av(n[1]), av(n[0]))), format="csr") - return self._aveE2CC + if(self.dim == 1): + return self.aveEx2CC + elif(self.dim == 2): + return 0.5*sp.hstack((self.aveEx2CC, self.aveEy2CC), format="csr") + elif(self.dim == 3): + return (1./3)*sp.hstack((self.aveEx2CC, self.aveEy2CC, self.aveEz2CC), format="csr") @property def aveE2CCV(self): "Construct the averaging operator on cell edges to cell centers." - if getattr(self, '_aveE2CCV', None) is None: + if(self.dim == 1): + return self.aveEx2CC + elif(self.dim == 2): + return sp.block_diag((self.aveEx2CC, self.aveEy2CC), format="csr") + elif(self.dim == 3): + return sp.block_diag((self.aveEx2CC, self.aveEy2CC, self.aveEz2CC), format="csr") + + @property + def aveEx2CC(self): + "Construct the averaging operator on cell edges in the x direction to cell centers." + if getattr(self, '_aveEx2CC', None) is None: # The number of cell centers in each direction n = self.vnC if(self.dim == 1): - raise Exception('Edge Averaging does not make sense in 1D: Use Identity?') + self._aveEx2CC = speye(n[0]) elif(self.dim == 2): - self._aveE2CCV = sp.block_diag((sp.kron(av(n[1]), speye(n[0])), - sp.kron(speye(n[1]), av(n[0]))), format="csr") + self._aveEx2CC = sp.kron(av(n[1]), speye(n[0])) elif(self.dim == 3): - self._aveE2CCV = sp.block_diag((kron3(av(n[2]), av(n[1]), speye(n[0])), - kron3(av(n[2]), speye(n[1]), av(n[0])), - kron3(speye(n[2]), av(n[1]), av(n[0]))), format="csr") - return self._aveE2CCV + self._aveEx2CC = kron3(av(n[2]), av(n[1]), speye(n[0])) + return self._aveEx2CC + + @property + def aveEy2CC(self): + "Construct the averaging operator on cell edges in the y direction to cell centers." + if self.dim < 2: return None + if getattr(self, '_aveEy2CC', None) is None: + # The number of cell centers in each direction + n = self.vnC + if(self.dim == 2): + self._aveEy2CC = sp.kron(speye(n[1]), av(n[0])) + elif(self.dim == 3): + self._aveEy2CC = kron3(av(n[2]), speye(n[1]), av(n[0])) + return self._aveEy2CC + + @property + def aveEz2CC(self): + "Construct the averaging operator on cell edges in the z direction to cell centers." + if self.dim < 3: return None + if getattr(self, '_aveEz2CC', None) is None: + # The number of cell centers in each direction + n = self.vnC + if(self.dim == 3): + self._aveEz2CC = kron3(speye(n[2]), av(n[1]), av(n[0])) + return self._aveEz2CC @property def aveN2CC(self): From b518b9b871f4d95251b9666c6d6f24f74358b38e Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Thu, 19 Mar 2015 16:22:55 -0700 Subject: [PATCH 3/3] updates to cylindrical to tensor mesh interpolation, tested. --- SimPEG/Mesh/CylMesh.py | 58 +++++++++++++++++++++++++++- SimPEG/Mesh/View.py | 23 +++++++++-- SimPEG/Tests/test_cylMesh.py | 74 ++++++++++++++++++++++++++++++++++++ 3 files changed, 151 insertions(+), 4 deletions(-) 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']