From 10ad9a762a8c097ed7d53037bc234db5fcf1cdcd Mon Sep 17 00:00:00 2001 From: Lindsey Date: Tue, 24 Feb 2015 13:49:33 -0800 Subject: [PATCH 1/7] test for Probing diagonal estimator. --- SimPEG/Tests/TestUtils.py | 2 +- SimPEG/Tests/test_utils.py | 13 +++++++++++++ 2 files changed, 14 insertions(+), 1 deletion(-) diff --git a/SimPEG/Tests/TestUtils.py b/SimPEG/Tests/TestUtils.py index 3cd40991..4292b0b4 100644 --- a/SimPEG/Tests/TestUtils.py +++ b/SimPEG/Tests/TestUtils.py @@ -1,7 +1,7 @@ import numpy as np import matplotlib.pyplot as plt from numpy.linalg import norm -from SimPEG.Utils import mkvc, sdiag +from SimPEG.Utils import mkvc, sdiag, diagEst from SimPEG import Utils from SimPEG.Mesh import TensorMesh, LogicallyRectMesh, CylMesh import numpy as np diff --git a/SimPEG/Tests/test_utils.py b/SimPEG/Tests/test_utils.py index 50cc0574..fe9e07c4 100644 --- a/SimPEG/Tests/test_utils.py +++ b/SimPEG/Tests/test_utils.py @@ -246,5 +246,18 @@ class TestSequenceFunctions(unittest.TestCase): self.assertTrue(np.all(true == listArray)) self.assertTrue(true.shape == listArray.shape) +class TestDiagEst(unittest.TestCase): + + def setUp(self): + self.n = 10 + self.A = np.random.rand(self.n,self.n) + self.Adiag = np.diagonal(self.A) + + def testOnes(self): + Adiagtest = diagEst(self.A,self.n,self.n) + r = np.abs(Adiagtest-self.Adiag) + self.assertTrue(r.dot(r) < TOL) + + if __name__ == '__main__': unittest.main() From 315bf7b61b359a19cba2ada7606d2b43e2b73955 Mon Sep 17 00:00:00 2001 From: Lindsey Date: Tue, 24 Feb 2015 15:34:41 -0800 Subject: [PATCH 2/7] added documentation to diagEst --- SimPEG/Utils/matutils.py | 31 +++++++++++++++++++++++++------ 1 file changed, 25 insertions(+), 6 deletions(-) diff --git a/SimPEG/Utils/matutils.py b/SimPEG/Utils/matutils.py index 57d165ad..e4ab5700 100644 --- a/SimPEG/Utils/matutils.py +++ b/SimPEG/Utils/matutils.py @@ -341,30 +341,49 @@ def invPropertyTensor(M, tensor, returnMatrix=False): return T -def diagEst(matFun, n, k=None, type='Probing'): - """ Based on Saad http://www-users.cs.umn.edu/~saad/PDF/umsi-2005-082.pdf, and http://www.cita.utoronto.ca/~niels/diagonal.pdf""" +def diagEst(matFun, n, k=None, approach='Probing'): + """ + Estimate the diagonal of a matrix, A. Note that the matrix may be a function which returns A times a vector. + + Three different approaches have been implemented, + 1. Probing : uses cyclic permutations of vectors with ones and zeros (default) + 2. Ones : random +/- 1 entries + 3. Random : random vectors + + :param lambda (numpy.array) matFun: matrix to estimate the diagonal of + :param int64 n: size of the vector that should be used to compute matFun(v) + :param int64 k: number of vectors to be used to estimate the diagonal + :param str approach: approach to be used for getting vectors + :rtype: numpy.array + :return: est_diag(A) + + Based on Saad http://www-users.cs.umn.edu/~saad/PDF/umsi-2005-082.pdf, and http://www.cita.utoronto.ca/~niels/diagonal.pdf + """ + + if type(matFun).__name__=='ndarray': + A = matFun + matFun = lambda v: A.dot(v) if k is None: k = np.floor(n/10.) - if type =='Ones': + if approach =='Ones': def getv(n,i=None): v = np.random.randn(n) v[v<0] = -1. v[v>=0] = 1. return v - elif type == 'Random': + elif approach == 'Random': def getv(n,i=None): return np.random.randn(n) - else: #if type == 'Probing': + else: #if approach == 'Probing': def getv(n,i): v = np.zeros(n) v[i:n:k] = 1. return v - #d = np.zeros(n) Mv = np.zeros(n) vv = np.zeros(n) From d22e72d0a7abcd0f009e862686721be1b91d5d27 Mon Sep 17 00:00:00 2001 From: Lindsey Date: Wed, 25 Feb 2015 16:09:02 -0800 Subject: [PATCH 3/7] If E0 is below the tolerance, the deriv test should pass --- SimPEG/Tests/TestUtils.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SimPEG/Tests/TestUtils.py b/SimPEG/Tests/TestUtils.py index 3cd40991..08f1fa60 100644 --- a/SimPEG/Tests/TestUtils.py +++ b/SimPEG/Tests/TestUtils.py @@ -267,7 +267,7 @@ def checkDerivative(fctn, x0, num=7, plotIt=True, dx=None, expectedOrder=2, tole # Ensure we are about precision order0 = order0[E0[1:] > eps] order1 = order1[E1[1:] > eps] - belowTol = order1.size == 0 and order0.size > 0 + belowTol = (order1.size == 0 and order0.size >= 0) # Make sure we get the correct order correctOrder = order1.size > 0 and np.mean(order1) > tolerance * expectedOrder From 90b3dfa61d33214ca27b2f123785f85d0eb76565 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Tue, 17 Mar 2015 14:07:52 -0700 Subject: [PATCH 4/7] Use meshTensor to create the time mesh in a TD problem. --- SimPEG/Problem.py | 14 +------------- 1 file changed, 1 insertion(+), 13 deletions(-) diff --git a/SimPEG/Problem.py b/SimPEG/Problem.py index dea83b1a..a4251a8f 100644 --- a/SimPEG/Problem.py +++ b/SimPEG/Problem.py @@ -437,19 +437,7 @@ class BaseTimeProblem(BaseProblem): del self.timeMesh return - if type(value) is not list: - raise Exception('timeSteps must be a np.ndarray or a list of scalars and tuples.') - - proposed = [] - for v in value: - if Utils.isScalar(v): - proposed += [float(v)] - elif type(v) is tuple and len(v) == 2: - proposed += [float(v[0])]*int(v[1]) - else: - raise Exception('timeSteps list must contain only scalars and len(2) tuples.') - - self._timeSteps = np.array(proposed) + self._timeSteps = Utils.meshTensor(value) del self.timeMesh @property From 16c269e01350b856e3e64e36cd200b7530167851 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Thu, 19 Mar 2015 12:48:50 -0700 Subject: [PATCH 5/7] 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 6/7] 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 7/7] 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']