From 9512b377c0cf8199052621a8a18b32aa639fb023 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Mon, 16 Jun 2014 12:17:01 -0600 Subject: [PATCH] InnerProducts working as an operator. Simplifications and generalizations in inner product code. --- SimPEG/Mesh/InnerProducts.py | 257 +++++++++++------------- SimPEG/Mesh/TensorMesh.py | 97 +++------ SimPEG/Tests/test_massMatricesDerivs.py | 117 ++++------- SimPEG/Tests/test_utils.py | 16 +- SimPEG/Utils/matutils.py | 57 ++++-- 5 files changed, 230 insertions(+), 314 deletions(-) diff --git a/SimPEG/Mesh/InnerProducts.py b/SimPEG/Mesh/InnerProducts.py index 8268b29a..026e406d 100644 --- a/SimPEG/Mesh/InnerProducts.py +++ b/SimPEG/Mesh/InnerProducts.py @@ -10,112 +10,43 @@ class InnerProducts(object): def __init__(self): raise Exception('InnerProducts is a base class providing inner product matrices for meshes and cannot run on its own. Inherit to your favorite Mesh class.') - def getFaceInnerProduct(self, prop=None, returnP=False, - invProp=False, invMat=False, doFast=True): + def getFaceInnerProduct(self, prop=None, invProp=False, invMat=False, doFast=True): """ :param numpy.array prop: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6)) - :param bool returnP: returns the projection matrices :param bool invProp: inverts the material property :param bool invMat: inverts the matrix :param bool doFast: do a faster implementation if available. :rtype: scipy.csr_matrix :return: M, the inner product matrix (nF, nF) """ - fast = None + return self._getInnerProduct('F', prop=prop, invProp=invProp, invMat=invMat, doFast=True) - if returnP is False and hasattr(self, '_fastFaceInnerProduct') and doFast: - fast = self._fastFaceInnerProduct(prop=prop, invProp=invProp, invMat=invMat) - - if fast is not None: - return fast - - if invProp: - prop = invPropertyTensor(self, prop) - - Mu = makePropertyTensor(self, prop) - - d = self.dim - # We will multiply by sqrt on each side to keep symmetry - V = sp.kron(sp.identity(d), sdiag(np.sqrt((2**(-d))*self.vol))) - - if d == 1: - fP = _getFacePx(self) - P000 = V*fP('fXm') - P100 = V*fP('fXp') - elif d == 2: - fP = _getFacePxx(self) - P000 = V*fP('fXm', 'fYm') - P100 = V*fP('fXp', 'fYm') - P010 = V*fP('fXm', 'fYp') - P110 = V*fP('fXp', 'fYp') - elif d == 3: - fP = _getFacePxxx(self) - P000 = V*fP('fXm', 'fYm', 'fZm') - P100 = V*fP('fXp', 'fYm', 'fZm') - P010 = V*fP('fXm', 'fYp', 'fZm') - P110 = V*fP('fXp', 'fYp', 'fZm') - P001 = V*fP('fXm', 'fYm', 'fZp') - P101 = V*fP('fXp', 'fYm', 'fZp') - P011 = V*fP('fXm', 'fYp', 'fZp') - P111 = V*fP('fXp', 'fYp', 'fZp') - - A = P000.T*Mu*P000 + P100.T*Mu*P100 - P = [P000, P100] - - if d > 1: - A = A + P010.T*Mu*P010 + P110.T*Mu*P110 - P += [P010, P110] - if d > 2: - A = A + P001.T*Mu*P001 + P101.T*Mu*P101 + P011.T*Mu*P011 + P111.T*Mu*P111 - P += [P001, P101, P011, P111] - - if invMat and tensorType(self, prop) < 3: - A = sdInv(A) - elif invMat and tensorType(self, prop) == 3: - raise Exception('Solver needed to invert A.') - - if returnP: - return A, P - else: - return A - - def getFaceInnerProductDeriv(self, prop=None, v=None, P=None, doFast=True): + def getEdgeInnerProduct(self, prop=None, invProp=False, invMat=False, doFast=True): """ :param numpy.array prop: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6)) - :param numpy.array v: vector to multiply (required in the general implementation) - :param list P: list of projection matrices - :param bool doFast: do a faster implementation if available. - :rtype: scipy.csr_matrix - :return: dMdm, the derivative of the inner product matrix (nF, nC*nA) - """ - fast = None - - if hasattr(self, '_fastFaceInnerProductDeriv') and doFast: - fast = self._fastFaceInnerProductDeriv(prop=prop, v=v) - - if fast is not None: - return fast - - if P is None: - M, P = self.getFaceInnerProduct(prop=prop, returnP=True) - - return self._getInnerProductDeriv(prop, v, P, self.nF) - - def getEdgeInnerProduct(self, prop=None, returnP=False, - invProp=False, invMat=False, doFast=True): - """ - :param numpy.array prop: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6)) - :param bool returnP: returns the projection matrices :param bool invProp: inverts the material property :param bool invMat: inverts the matrix :param bool doFast: do a faster implementation if available. :rtype: scipy.csr_matrix :return: M, the inner product matrix (nE, nE) """ + return self._getInnerProduct('E', prop=prop, invProp=invProp, invMat=invMat) + + def _getInnerProduct(self, projType, prop=None, invProp=False, invMat=False, doFast=True): + """ + :param str projType: 'F' for faces 'E' for edges + :param numpy.array prop: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6)) + :param bool invProp: inverts the material property + :param bool invMat: inverts the matrix + :param bool doFast: do a faster implementation if available. + :rtype: scipy.csr_matrix + :return: M, the inner product matrix (nE, nE) + """ + assert projType in ['F', 'E'], "projType must be 'F' for faces or 'E' for edges" fast = None - if returnP is False and hasattr(self, '_fastEdgeInnerProduct') and doFast: - fast = self._fastEdgeInnerProduct(prop=prop, invProp=invProp, invMat=invMat) + if hasattr(self, '_fastInnerProduct') and doFast: + fast = self._fastInnerProduct(projType, prop=prop, invProp=invProp, invMat=invMat) if fast is not None: return fast @@ -123,72 +54,122 @@ class InnerProducts(object): if invProp: prop = invPropertyTensor(self, prop) + tensorType = TensorType(self, prop) Mu = makePropertyTensor(self, prop) + Ps = self._getInnerProductProjectionMatrices(projType, tensorType) + + A = np.sum([P.T * Mu * P for P in Ps]) + + if invMat and tensorType < 3: + A = sdInv(A) + elif invMat and tensorType == 3: + raise Exception('Solver needed to invert A.') + + return A + + def _getInnerProductProjectionMatrices(self, projType, tensorType): + """ + :param str projType: 'F' for faces 'E' for edges + :param TensorType tensorType: type of the tensor: TensorType(mesh, sigma) + """ + assert isinstance(tensorType, TensorType), 'tensorType must be an instance of TensorType.' + assert projType in ['F', 'E'], "projType must be 'F' for faces or 'E' for edges" + d = self.dim # We will multiply by sqrt on each side to keep symmetry V = sp.kron(sp.identity(d), sdiag(np.sqrt((2**(-d))*self.vol))) - if d == 1: - raise NotImplementedError('getEdgeInnerProduct not implemented for 1D') - elif d == 2: - eP = _getEdgePxx(self) - P000 = V*eP('eX0', 'eY0') - P100 = V*eP('eX0', 'eY1') - P010 = V*eP('eX1', 'eY0') - P110 = V*eP('eX1', 'eY1') - elif d == 3: - eP = _getEdgePxxx(self) - P000 = V*eP('eX0', 'eY0', 'eZ0') - P100 = V*eP('eX0', 'eY1', 'eZ1') - P010 = V*eP('eX1', 'eY0', 'eZ2') - P110 = V*eP('eX1', 'eY1', 'eZ3') - P001 = V*eP('eX2', 'eY2', 'eZ0') - P101 = V*eP('eX2', 'eY3', 'eZ1') - P011 = V*eP('eX3', 'eY2', 'eZ2') - P111 = V*eP('eX3', 'eY3', 'eZ3') + nodes = ['000', '100', '010', '110', '001', '101', '011', '111'][:2**d] - Mu = makePropertyTensor(self, prop) - A = P000.T*Mu*P000 + P100.T*Mu*P100 + P010.T*Mu*P010 + P110.T*Mu*P110 - P = [P000, P100, P010, P110] - if d == 3: - A = A + P001.T*Mu*P001 + P101.T*Mu*P101 + P011.T*Mu*P011 + P111.T*Mu*P111 - P += [P001, P101, P011, P111] + if projType == 'F': + locs = { + '000': [('fXm',), ('fXm', 'fYm'), ('fXm', 'fYm', 'fZm')], + '100': [('fXp',), ('fXp', 'fYm'), ('fXp', 'fYm', 'fZm')], + '010': [ None , ('fXm', 'fYp'), ('fXm', 'fYp', 'fZm')], + '110': [ None , ('fXp', 'fYp'), ('fXp', 'fYp', 'fZm')], + '001': [ None , None , ('fXm', 'fYm', 'fZp')], + '101': [ None , None , ('fXp', 'fYm', 'fZp')], + '011': [ None , None , ('fXm', 'fYp', 'fZp')], + '111': [ None , None , ('fXp', 'fYp', 'fZp')] + } + if d == 1: + proj = _getFacePx(self) + elif d == 2: + proj = _getFacePxx(self) + elif d == 3: + proj = _getFacePxxx(self) - if invMat and tensorType(self, prop) < 3: - A = sdInv(A) - elif invMat and tensorType(self, prop) == 3: - raise Exception('Solver needed to invert A.') + elif projType == 'E': + locs = { + '000': [ None , ('eX0', 'eY0'), ('eX0', 'eY0', 'eZ0')], + '100': [ None , ('eX0', 'eY1'), ('eX0', 'eY1', 'eZ1')], + '010': [ None , ('eX1', 'eY0'), ('eX1', 'eY0', 'eZ2')], + '110': [ None , ('eX1', 'eY1'), ('eX1', 'eY1', 'eZ3')], + '001': [ None , None , ('eX2', 'eY2', 'eZ0')], + '101': [ None , None , ('eX2', 'eY3', 'eZ1')], + '011': [ None , None , ('eX3', 'eY2', 'eZ2')], + '111': [ None , None , ('eX3', 'eY3', 'eZ3')] + } + if d == 1: + raise NotImplementedError('getEdgeInnerProduct not implemented for 1D') + elif d == 2: + proj = _getEdgePxx(self) + elif d == 3: + proj = _getEdgePxxx(self) - if returnP: - return A, P - else: - return A + return [V*proj(*locs[node][d-1]) for node in nodes] - def getEdgeInnerProductDeriv(self, prop=None, v=None, P=None, doFast=True): + + def getFaceInnerProductDeriv(self, tensorType, P=None, doFast=True): """ - :param numpy.array prop: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6)) - :param numpy.array v: vector to multiply (required in the general implementation) + :param TensorType tensorType: type of the tensor: TensorType(mesh, sigma) :param list P: list of projection matrices :param bool doFast: do a faster implementation if available. :rtype: scipy.csr_matrix - :return: dMdm, the derivative of the inner product matrix (nE, nC*nA) + :return: dMdm, the derivative of the inner product matrix (nF, nC*nA) """ - + assert isinstance(tensorType, TensorType), 'tensorType must be an instance of TensorType.' fast = None - if hasattr(self, '_fastEdgeInnerProductDeriv') and doFast: - fast = self._fastEdgeInnerProductDeriv(prop=prop, v=v) + if hasattr(self, '_fastInnerProductDeriv') and doFast: + fast = self._fastInnerProductDeriv('F', tensorType) if fast is not None: return fast if P is None: - M, P = self.getEdgeInnerProduct(prop=prop, returnP=True) + P = self._getInnerProductProjectionMatrices('F', tensorType=tensorType) - return self._getInnerProductDeriv(prop, v, P, self.nE) + def innerProductDeriv(v): + return self._getInnerProductDeriv(tensorType, P, self.nF, v) + return DerivOperator(innerProductDeriv) - def _getInnerProductDeriv(self, prop, v, P, n): + + def getEdgeInnerProductDeriv(self, tensorType, P=None, doFast=True): + """ + :param TensorType tensorType: type of the tensor: TensorType(mesh, sigma) + :param list P: list of projection matrices + :param bool doFast: do a faster implementation if available. + :rtype: scipy.csr_matrix + :return: dMdm, the derivative of the inner product matrix (nE, nC*nA) + """ + assert isinstance(tensorType, TensorType), 'tensorType must be an instance of TensorType.' + fast = None + + if hasattr(self, '_fastInnerProductDeriv') and doFast: + fast = self._fastInnerProductDeriv('E', tensorType) + + if fast is not None: + return fast + + if P is None: + P = self._getInnerProductProjectionMatrices('E', tensorType=tensorType) + def innerProductDeriv(v): + return self._getInnerProductDeriv(tensorType, P, self.nE, v) + return DerivOperator(innerProductDeriv) + + def _getInnerProductDeriv(self, tensorType, P, n, v): """ :param numpy.array prop: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6)) :param numpy.array v: vector to multiply (required in the general implementation) @@ -197,7 +178,7 @@ class InnerProducts(object): :rtype: scipy.csr_matrix :return: dMdm, the derivative of the inner product matrix (n, nC*nA) """ - if prop is None: + if tensorType == -1: return None if v is None: @@ -206,24 +187,24 @@ class InnerProducts(object): d = self.dim Z = spzeros(self.nC, self.nC) - if isScalar(prop): + if tensorType == 0: dMdm = spzeros(n, 1) for i, p in enumerate(P): dMdm = dMdm + sp.csr_matrix((p.T * (p * v), (range(n), np.zeros(n))), shape=(n,1)) if d == 1: - if prop.size == self.nC: + if tensorType == 1: dMdm = spzeros(n, self.nC) for i, p in enumerate(P): dMdm = dMdm + p.T * sdiag( p * v ) elif d == 2: - if prop.size == self.nC: + if tensorType == 1: dMdm = spzeros(n, self.nC) for i, p in enumerate(P): Y = p * v y1 = Y[:self.nC] y2 = Y[self.nC:] dMdm = dMdm + p.T * sp.vstack((sdiag( y1 ), sdiag( y2 ))) - elif prop.size == self.nC*2: + elif tensorType == 2: dMdms = [spzeros(n, self.nC) for _ in range(2)] for i, p in enumerate(P): Y = p * v @@ -232,7 +213,7 @@ class InnerProducts(object): dMdms[0] = dMdms[0] + p.T * sp.vstack(( sdiag( y1 ), Z)) dMdms[1] = dMdms[1] + p.T * sp.vstack(( Z, sdiag( y2 ))) dMdm = sp.hstack(dMdms) - elif prop.size == self.nC*3: + elif tensorType == 3: dMdms = [spzeros(n, self.nC) for _ in range(3)] for i, p in enumerate(P): Y = p * v @@ -243,7 +224,7 @@ class InnerProducts(object): dMdms[2] = dMdms[2] + p.T * sp.vstack(( sdiag( y2 ), sdiag( y1 ))) dMdm = sp.hstack(dMdms) elif d == 3: - if prop.size == self.nC: + if tensorType == 1: dMdm = spzeros(n, self.nC) for i, p in enumerate(P): Y = p * v @@ -251,7 +232,7 @@ class InnerProducts(object): y2 = Y[self.nC:self.nC*2] y3 = Y[self.nC*2:] dMdm = dMdm + p.T * sp.vstack((sdiag( y1 ), sdiag( y2 ), sdiag( y3 ))) - elif prop.size == self.nC*3: + elif tensorType == 2: dMdms = [spzeros(n, self.nC) for _ in range(3)] for i, p in enumerate(P): Y = p * v @@ -262,7 +243,7 @@ class InnerProducts(object): dMdms[1] = dMdms[1] + p.T * sp.vstack(( Z, sdiag( y2 ), Z)) dMdms[2] = dMdms[2] + p.T * sp.vstack(( Z, Z, sdiag( y3 ))) dMdm = sp.hstack(dMdms) - elif prop.size == self.nC*6: + elif tensorType == 3: dMdms = [spzeros(n, self.nC) for _ in range(6)] for i, p in enumerate(P): Y = p * v diff --git a/SimPEG/Mesh/TensorMesh.py b/SimPEG/Mesh/TensorMesh.py index ba5bd258..faf45cf2 100644 --- a/SimPEG/Mesh/TensorMesh.py +++ b/SimPEG/Mesh/TensorMesh.py @@ -260,49 +260,21 @@ class BaseTensorMesh(BaseRectangularMesh): return Q.tocsr() - def _fastFaceInnerProduct(self, prop=None, invProp=False, invMat=False): + def _fastInnerProduct(self, projType, prop=None, invProp=False, invMat=False): """ Fast version of getFaceInnerProduct. This does not handle the case of a full tensor prop. :param numpy.array prop: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6)) + :param str projType: 'E' or 'F' :param bool returnP: returns the projection matrices :param bool invProp: inverts the material property :param bool invMat: inverts the matrix :rtype: scipy.csr_matrix :return: M, the inner product matrix (nF, nF) """ - return self._fastInnerProduct('F', prop=prop, invProp=invProp, invMat=invMat) + assert projType in ['F', 'E'], "projType must be 'F' for faces or 'E' for edges" - - def _fastEdgeInnerProduct(self, prop=None, invProp=False, invMat=False): - """ - Fast version of getEdgeInnerProduct. - This does not handle the case of a full tensor prop. - - :param numpy.array prop: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6)) - :param bool returnP: returns the projection matrices - :param bool invProp: inverts the material property - :param bool invMat: inverts the matrix - :rtype: scipy.csr_matrix - :return: M, the inner product matrix (nE, nE) - """ - return self._fastInnerProduct('E', prop=prop, invProp=invProp, invMat=invMat) - - - def _fastInnerProduct(self, AvType, prop=None, invProp=False, invMat=False): - """ - Fast version of getFaceInnerProduct. - This does not handle the case of a full tensor prop. - - :param numpy.array prop: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6)) - :param str AvType: 'E' or 'F' - :param bool returnP: returns the projection matrices - :param bool invProp: inverts the material property - :param bool invMat: inverts the matrix - :rtype: scipy.csr_matrix - :return: M, the inner product matrix (nF, nF) - """ if prop is None: prop = np.ones(self.nC) @@ -313,11 +285,11 @@ class BaseTensorMesh(BaseRectangularMesh): prop = prop*np.ones(self.nC) if prop.size == self.nC: - Av = getattr(self, 'ave'+AvType+'2CC') + Av = getattr(self, 'ave'+projType+'2CC') Vprop = self.vol * Utils.mkvc(prop) M = self.dim * Utils.sdiag(Av.T * Vprop) elif prop.size == self.nC*self.dim: - Av = getattr(self, 'ave'+AvType+'2CCV') + Av = getattr(self, 'ave'+projType+'2CCV') V = sp.kron(sp.identity(self.dim), Utils.sdiag(self.vol)) M = Utils.sdiag(Av.T * V * Utils.mkvc(prop)) else: @@ -328,55 +300,40 @@ class BaseTensorMesh(BaseRectangularMesh): else: return M - def _fastFaceInnerProductDeriv(self, prop=None, v=None): + def _fastInnerProductDeriv(self, projType, tensorType): """ - :param numpy.array prop: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6)) + :param str projType: 'E' or 'F' + :param TensorType tensorType: type of the tensor :rtype: scipy.csr_matrix :return: M, the inner product matrix (nF, nF) """ - return self._fastInnerProductDeriv('F', prop=prop, v=v) - - - def _fastEdgeInnerProductDeriv(self, prop=None, v=None): - """ - :param numpy.array prop: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6)) - :rtype: scipy.csr_matrix - :return: M, the inner product matrix (nE, nE) - """ - return self._fastInnerProductDeriv('E', prop=prop, v=v) - - - def _fastInnerProductDeriv(self, AvType, prop=None, v=None): - """ - :param str AvType: 'E' or 'F' - :param numpy.array prop: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6)) - :rtype: scipy.csr_matrix - :return: M, the inner product matrix (nF, nF) - """ - if prop is None: + assert projType in ['F', 'E'], "projType must be 'F' for faces or 'E' for edges" + if tensorType == -1: return None - if Utils.isScalar(prop): - Av = getattr(self, 'ave'+AvType+'2CC') + if tensorType == 0: + Av = getattr(self, 'ave'+projType+'2CC') V = Utils.sdiag(self.vol) ones = sp.csr_matrix((np.ones(self.nC), (range(self.nC), np.zeros(self.nC))), shape=(self.nC,1)) - if v is None: - return self.dim * Av.T * V * ones - return Utils.sdiag(v) * self.dim * Av.T * V * ones + # if v is None: + # return self.dim * Av.T * V * ones + def scalarInnerProductDeriv(v): + return Utils.sdiag(v) * self.dim * Av.T * V * ones + return Utils.DerivOperator(scalarInnerProductDeriv) - if prop.size == self.nC: - Av = getattr(self, 'ave'+AvType+'2CC') + if tensorType == 1: + Av = getattr(self, 'ave'+projType+'2CC') V = Utils.sdiag(self.vol) - if v is None: - return self.dim * Av.T * V - return Utils.sdiag(v) * self.dim * Av.T * V + def isotropicInnerProductDeriv(v): + return Utils.sdiag(v) * self.dim * Av.T * V + return Utils.DerivOperator(isotropicInnerProductDeriv) - if prop.size == self.nC*self.dim: # anisotropic - Av = getattr(self, 'ave'+AvType+'2CCV') + if tensorType == 2: # anisotropic + Av = getattr(self, 'ave'+projType+'2CCV') V = sp.kron(sp.identity(self.dim), Utils.sdiag(self.vol)) - if v is None: - return Av.T * V - return Utils.sdiag(v) * Av.T * V + def anisotropicInnerProductDeriv(v): + return Utils.sdiag(v) * Av.T * V + return Utils.DerivOperator(anisotropicInnerProductDeriv) diff --git a/SimPEG/Tests/test_massMatricesDerivs.py b/SimPEG/Tests/test_massMatricesDerivs.py index 1da88173..b574c8f1 100644 --- a/SimPEG/Tests/test_massMatricesDerivs.py +++ b/SimPEG/Tests/test_massMatricesDerivs.py @@ -6,130 +6,93 @@ from TestUtils import checkDerivative class TestInnerProductsDerivs(unittest.TestCase): - def doTestFace(self, h, rep, vec, fast): + def doTestFace(self, h, rep, fast): mesh = Mesh.TensorMesh(h) v = np.random.rand(mesh.nF) + sig = np.random.rand(1) if rep is 0 else np.random.rand(mesh.nC*rep) + Md = mesh.getFaceInnerProductDeriv(Utils.TensorType(mesh, sig), doFast=fast) def fun(sig): M = mesh.getFaceInnerProduct(sig) - if vec: - Md = mesh.getFaceInnerProductDeriv(sig, v=v, doFast=fast) - return M*v, Md - Md = mesh.getFaceInnerProductDeriv(sig, doFast=fast) - return M*v, Utils.sdiag(v)*Md - sig = np.random.rand(1) if rep is 0 else np.random.rand(mesh.nC*rep) + return M*v, Md*v return checkDerivative(fun, sig, num=5, plotIt=False) - def doTestEdge(self, h, rep, vec, fast): + def doTestEdge(self, h, rep, fast): mesh = Mesh.TensorMesh(h) v = np.random.rand(mesh.nE) + sig = np.random.rand(1) if rep is 0 else np.random.rand(mesh.nC*rep) + Md = mesh.getEdgeInnerProductDeriv(Utils.TensorType(mesh, sig), doFast=fast) def fun(sig): M = mesh.getEdgeInnerProduct(sig) - if vec: - Md = mesh.getEdgeInnerProductDeriv(sig, v=v, doFast=fast) - return M*v, Md - Md = mesh.getEdgeInnerProductDeriv(sig, doFast=fast) - return M*v, Utils.sdiag(v)*Md - sig = np.random.rand(1) if rep is 0 else np.random.rand(mesh.nC*rep) + return M*v, Md*v return checkDerivative(fun, sig, num=5, plotIt=False) def test_FaceIP_1D_float(self): - self.assertTrue(self.doTestFace([10],0,True, False)) + self.assertTrue(self.doTestFace([10],0, False)) def test_FaceIP_2D_float(self): - self.assertTrue(self.doTestFace([10, 4],0,True, False)) + self.assertTrue(self.doTestFace([10, 4],0, False)) def test_FaceIP_3D_float(self): - self.assertTrue(self.doTestFace([10, 4, 5],0,True, False)) + self.assertTrue(self.doTestFace([10, 4, 5],0, False)) def test_FaceIP_1D_isotropic(self): - self.assertTrue(self.doTestFace([10],1,True, False)) + self.assertTrue(self.doTestFace([10],1, False)) def test_FaceIP_2D_isotropic(self): - self.assertTrue(self.doTestFace([10, 4],1,True, False)) + self.assertTrue(self.doTestFace([10, 4],1, False)) def test_FaceIP_3D_isotropic(self): - self.assertTrue(self.doTestFace([10, 4, 5],1,True, False)) + self.assertTrue(self.doTestFace([10, 4, 5],1, False)) def test_FaceIP_2D_anisotropic(self): - self.assertTrue(self.doTestFace([10, 4],2,True, False)) + self.assertTrue(self.doTestFace([10, 4],2, False)) def test_FaceIP_3D_anisotropic(self): - self.assertTrue(self.doTestFace([10, 4, 5],3,True, False)) + self.assertTrue(self.doTestFace([10, 4, 5],3, False)) def test_FaceIP_2D_tensor(self): - self.assertTrue(self.doTestFace([10, 4],3,True, False)) + self.assertTrue(self.doTestFace([10, 4],3, False)) def test_FaceIP_3D_tensor(self): - self.assertTrue(self.doTestFace([10, 4, 5],6,True, False)) + self.assertTrue(self.doTestFace([10, 4, 5],6, False)) def test_FaceIP_1D_float_fast(self): - self.assertTrue(self.doTestFace([10],0, False, True)) + self.assertTrue(self.doTestFace([10],0, True)) def test_FaceIP_2D_float_fast(self): - self.assertTrue(self.doTestFace([10, 4],0, False, True)) + self.assertTrue(self.doTestFace([10, 4],0, True)) def test_FaceIP_3D_float_fast(self): - self.assertTrue(self.doTestFace([10, 4, 5],0, False, True)) + self.assertTrue(self.doTestFace([10, 4, 5],0, True)) def test_FaceIP_1D_isotropic_fast(self): - self.assertTrue(self.doTestFace([10],1, False, True)) + self.assertTrue(self.doTestFace([10],1, True)) def test_FaceIP_2D_isotropic_fast(self): - self.assertTrue(self.doTestFace([10, 4],1, False, True)) + self.assertTrue(self.doTestFace([10, 4],1, True)) def test_FaceIP_3D_isotropic_fast(self): - self.assertTrue(self.doTestFace([10, 4, 5],1, False, True)) + self.assertTrue(self.doTestFace([10, 4, 5],1, True)) def test_FaceIP_2D_anisotropic_fast(self): - self.assertTrue(self.doTestFace([10, 4],2, False, True)) + self.assertTrue(self.doTestFace([10, 4],2, True)) def test_FaceIP_3D_anisotropic_fast(self): - self.assertTrue(self.doTestFace([10, 4, 5],3, False, True)) - - def test_FaceIP_1D_float_fast_vec(self): - self.assertTrue(self.doTestFace([10],0, True, True)) - def test_FaceIP_2D_float_fast_vec(self): - self.assertTrue(self.doTestFace([10, 4],0, True, True)) - def test_FaceIP_3D_float_fast_vec(self): - self.assertTrue(self.doTestFace([10, 4, 5],0, True, True)) - def test_FaceIP_1D_isotropic_fast_vec(self): - self.assertTrue(self.doTestFace([10],1, True, True)) - def test_FaceIP_2D_isotropic_fast_vec(self): - self.assertTrue(self.doTestFace([10, 4],1, True, True)) - def test_FaceIP_3D_isotropic_fast_vec(self): - self.assertTrue(self.doTestFace([10, 4, 5],1, True, True)) - def test_FaceIP_2D_anisotropic_fast_vec(self): - self.assertTrue(self.doTestFace([10, 4],2, True, True)) - def test_FaceIP_3D_anisotropic_fast_vec(self): - self.assertTrue(self.doTestFace([10, 4, 5],3, True, True)) + self.assertTrue(self.doTestFace([10, 4, 5],3, True)) def test_EdgeIP_2D_float(self): - self.assertTrue(self.doTestEdge([10, 4],0,True, False)) + self.assertTrue(self.doTestEdge([10, 4],0, False)) def test_EdgeIP_3D_float(self): - self.assertTrue(self.doTestEdge([10, 4, 5],0,True, False)) + self.assertTrue(self.doTestEdge([10, 4, 5],0, False)) def test_EdgeIP_2D_isotropic(self): - self.assertTrue(self.doTestEdge([10, 4],1,True, False)) + self.assertTrue(self.doTestEdge([10, 4],1, False)) def test_EdgeIP_3D_isotropic(self): - self.assertTrue(self.doTestEdge([10, 4, 5],1,True, False)) + self.assertTrue(self.doTestEdge([10, 4, 5],1, False)) def test_EdgeIP_2D_anisotropic(self): - self.assertTrue(self.doTestEdge([10, 4],2,True, False)) + self.assertTrue(self.doTestEdge([10, 4],2, False)) def test_EdgeIP_3D_anisotropic(self): - self.assertTrue(self.doTestEdge([10, 4, 5],3,True, False)) + self.assertTrue(self.doTestEdge([10, 4, 5],3, False)) def test_EdgeIP_2D_tensor(self): - self.assertTrue(self.doTestEdge([10, 4],3,True, False)) + self.assertTrue(self.doTestEdge([10, 4],3, False)) def test_EdgeIP_3D_tensor(self): - self.assertTrue(self.doTestEdge([10, 4, 5],6,True, False)) + self.assertTrue(self.doTestEdge([10, 4, 5],6, False)) def test_EdgeIP_2D_float_fast(self): - self.assertTrue(self.doTestEdge([10, 4],0, False, True)) + self.assertTrue(self.doTestEdge([10, 4],0, True)) def test_EdgeIP_3D_float_fast(self): - self.assertTrue(self.doTestEdge([10, 4, 5],0, False, True)) + self.assertTrue(self.doTestEdge([10, 4, 5],0, True)) def test_EdgeIP_2D_isotropic_fast(self): - self.assertTrue(self.doTestEdge([10, 4],1, False, True)) + self.assertTrue(self.doTestEdge([10, 4],1, True)) def test_EdgeIP_3D_isotropic_fast(self): - self.assertTrue(self.doTestEdge([10, 4, 5],1, False, True)) + self.assertTrue(self.doTestEdge([10, 4, 5],1, True)) def test_EdgeIP_2D_anisotropic_fast(self): - self.assertTrue(self.doTestEdge([10, 4],2, False, True)) + self.assertTrue(self.doTestEdge([10, 4],2, True)) def test_EdgeIP_3D_anisotropic_fast(self): - self.assertTrue(self.doTestEdge([10, 4, 5],3, False, True)) - - def test_EdgeIP_2D_float_fast_vec(self): - self.assertTrue(self.doTestEdge([10, 4],0, True, True)) - def test_EdgeIP_3D_float_fast_vec(self): - self.assertTrue(self.doTestEdge([10, 4, 5],0, True, True)) - def test_EdgeIP_2D_isotropic_fast_vec(self): - self.assertTrue(self.doTestEdge([10, 4],1, True, True)) - def test_EdgeIP_3D_isotropic_fast_vec(self): - self.assertTrue(self.doTestEdge([10, 4, 5],1, True, True)) - def test_EdgeIP_2D_anisotropic_fast_vec(self): - self.assertTrue(self.doTestEdge([10, 4],2, True, True)) - def test_EdgeIP_3D_anisotropic_fast_vec(self): - self.assertTrue(self.doTestEdge([10, 4, 5],3, True, True)) - + self.assertTrue(self.doTestEdge([10, 4, 5],3, True)) if __name__ == '__main__': diff --git a/SimPEG/Tests/test_utils.py b/SimPEG/Tests/test_utils.py index 6318527a..50cc0574 100644 --- a/SimPEG/Tests/test_utils.py +++ b/SimPEG/Tests/test_utils.py @@ -160,7 +160,7 @@ class TestSequenceFunctions(unittest.TestCase): Z = B2*A - sp.identity(M.nC*2) self.assertTrue(np.linalg.norm(Z.todense().ravel(), 2) < TOL) - def test_tensorType2D(self): + def test_TensorType2D(self): M = Mesh.TensorMesh([6, 6]) a1 = np.random.rand(M.nC) a2 = np.random.rand(M.nC) @@ -170,12 +170,12 @@ class TestSequenceFunctions(unittest.TestCase): prop3 = np.c_[a1, a2, a3] for ii, prop in enumerate([4, prop1, prop2, prop3]): - self.assertTrue(tensorType(M, prop) == ii) + self.assertTrue(TensorType(M, prop) == ii) - self.assertRaises(Exception, tensorType, M, np.c_[a1, a2, a3, a3]) - self.assertTrue(tensorType(M, None) == -1) + self.assertRaises(Exception, TensorType, M, np.c_[a1, a2, a3, a3]) + self.assertTrue(TensorType(M, None) == -1) - def test_tensorType3D(self): + def test_TensorType3D(self): M = Mesh.TensorMesh([6, 6, 7]) a1 = np.random.rand(M.nC) a2 = np.random.rand(M.nC) @@ -188,10 +188,10 @@ class TestSequenceFunctions(unittest.TestCase): prop3 = np.c_[a1, a2, a3, a4, a5, a6] for ii, prop in enumerate([4, prop1, prop2, prop3]): - self.assertTrue(tensorType(M, prop) == ii) + self.assertTrue(TensorType(M, prop) == ii) - self.assertRaises(Exception, tensorType, M, np.c_[a1, a2, a3, a3]) - self.assertTrue(tensorType(M, None) == -1) + self.assertRaises(Exception, TensorType, M, np.c_[a1, a2, a3, a3]) + self.assertTrue(TensorType(M, None) == -1) def test_invPropertyTensor3D(self): diff --git a/SimPEG/Utils/matutils.py b/SimPEG/Utils/matutils.py index a93e6f07..1a4c4a0e 100644 --- a/SimPEG/Utils/matutils.py +++ b/SimPEG/Utils/matutils.py @@ -251,25 +251,34 @@ def inv2X2BlockDiagonal(a11, a12, a21, a22, returnMatrix=True): return sp.vstack((sp.hstack((sdiag(b11), sdiag(b12))), sp.hstack((sdiag(b21), sdiag(b22))))) -def tensorType(M, tensor): - if tensor is None: # default is ones - return -1 - - if isScalar(tensor): - return 0 - - if tensor.size == M.nC: - return 1 - - if ((M.dim == 2 and tensor.size == M.nC*2) or - (M.dim == 3 and tensor.size == M.nC*3)): - return 2 - - if ((M.dim == 2 and tensor.size == M.nC*3) or - (M.dim == 3 and tensor.size == M.nC*6)): - return 3 - - raise Exception('Unexpected shape of tensor') +class TensorType(object): + def __init__(self, M, tensor): + if tensor is None: # default is ones + self._tt = -1 + self._tts = 'none' + elif isScalar(tensor): + self._tt = 0 + self._tts = 'scalar' + elif tensor.size == M.nC: + self._tt = 1 + self._tts = 'isotropic' + elif ((M.dim == 2 and tensor.size == M.nC*2) or + (M.dim == 3 and tensor.size == M.nC*3)): + self._tt = 2 + self._tts = 'anisotropic' + elif ((M.dim == 2 and tensor.size == M.nC*3) or + (M.dim == 3 and tensor.size == M.nC*6)): + self._tt = 3 + self._tts = 'tensor' + else: + raise Exception('Unexpected shape of tensor') + def __str__(self): + return 'TensorType[%i]: %s' % (self._tt, self._tts) + def __eq__(self, v): return self._tt == v + def __le__(self, v): return self._tt <= v + def __ge__(self, v): return self._tt >= v + def __lt__(self, v): return self._tt < v + def __gt__(self, v): return self._tt > v def makePropertyTensor(M, tensor): if tensor is None: # default is ones @@ -278,7 +287,7 @@ def makePropertyTensor(M, tensor): if isScalar(tensor): tensor = tensor * np.ones(M.nC) - propType = tensorType(M, tensor) + propType = TensorType(M, tensor) if propType == 1: # Isotropic! Sigma = sp.kron(sp.identity(M.dim), sdiag(mkvc(tensor))) elif propType == 2: # Diagonal tensor @@ -302,7 +311,7 @@ def makePropertyTensor(M, tensor): def invPropertyTensor(M, tensor, returnMatrix=False): - propType = tensorType(M, tensor) + propType = TensorType(M, tensor) if isScalar(tensor): T = 1./tensor @@ -341,3 +350,9 @@ class SimPEGLinearOperator(LinearOperator): def T(self): return self.__class__((self.shape[1],self.shape[0]),self.rmatvec,rmatvec=self.matvec,matmat=self.matmat) + +class DerivOperator(object): + def __init__(self, f): + self.f = f + def __mul__(self, v): + return self.f(v)