diff --git a/SimPEG/Mesh/InnerProducts.py b/SimPEG/Mesh/InnerProducts.py index 99490674..22c6a521 100644 --- a/SimPEG/Mesh/InnerProducts.py +++ b/SimPEG/Mesh/InnerProducts.py @@ -1,5 +1,5 @@ from scipy import sparse as sp -from SimPEG.Utils import sub2ind, ndgrid, mkvc, getSubArray, sdiag, inv3X3BlockDiagonal, inv2X2BlockDiagonal, makePropertyTensor, invPropertyTensor, spzeros +from SimPEG.Utils import sub2ind, ndgrid, mkvc, getSubArray, sdiag, inv3X3BlockDiagonal, inv2X2BlockDiagonal, makePropertyTensor, invPropertyTensor, spzeros, isScalar import numpy as np @@ -183,25 +183,33 @@ class InnerProducts(object): :rtype: scipy.csr_matrix :return: dMdm, the derivative of the inner product matrix (n, nC*nA) """ + if materialProperty is None: + return None + if v is None: raise Exception('v must be supplied for this implementation.') d = self.dim Z = spzeros(self.nC, self.nC) - if d == 1: - dMdm = spzeros(n, self.nC) + if isScalar(materialProperty): + dMdm = spzeros(n, 1) for i, p in enumerate(P): - dMdm = dMdm + p.T * sdiag( p * v ) + dMdm = dMdm + sp.csr_matrix((p.T * (p * v), (range(n), np.zeros(n))), shape=(n,1)) + if d == 1: + if materialProperty.size == self.nC: + dMdm = spzeros(n, self.nC) + for i, p in enumerate(P): + dMdm = dMdm + p.T * sdiag( p * v ) elif d == 2: - if materialProperty is None or materialProperty.size == self.nC: + if materialProperty.size == self.nC: 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 ))) - if materialProperty.size == self.nC*2: + elif materialProperty.size == self.nC*2: dMdms = [spzeros(n, self.nC) for _ in range(2)] for i, p in enumerate(P): Y = p * v @@ -210,7 +218,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) - if materialProperty.size == self.nC*3: + elif materialProperty.size == self.nC*3: dMdms = [spzeros(n, self.nC) for _ in range(3)] for i, p in enumerate(P): Y = p * v @@ -221,7 +229,7 @@ class InnerProducts(object): dMdms[2] = dMdms[2] + p.T * sp.vstack(( sdiag( y2 ), sdiag( y1 ))) dMdm = sp.hstack(dMdms) elif d == 3: - if materialProperty is None or materialProperty.size == self.nC: + if materialProperty.size == self.nC: dMdm = spzeros(n, self.nC) for i, p in enumerate(P): Y = p * v @@ -229,7 +237,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 ))) - if materialProperty.size == self.nC*3: + elif materialProperty.size == self.nC*3: dMdms = [spzeros(n, self.nC) for _ in range(3)] for i, p in enumerate(P): Y = p * v @@ -240,7 +248,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) - if materialProperty.size == self.nC*6: + elif materialProperty.size == self.nC*6: 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 df049970..1bc7cadb 100644 --- a/SimPEG/Mesh/TensorMesh.py +++ b/SimPEG/Mesh/TensorMesh.py @@ -534,18 +534,18 @@ class TensorMesh(BaseRectangularMesh, TensorView, DiffOperators, InnerProducts): """ if materialProperty is None: materialProperty = np.ones(self.nC) - elif type(materialProperty) in [float, int, long]: - materialProperty = materialProperty * np.ones(M.nC) + + if invertProperty: + materialProperty = 1./materialProperty + + if Utils.isScalar(materialProperty): + materialProperty = materialProperty*np.ones(self.nC) if materialProperty.size == self.nC: - if invertProperty: - materialProperty = 1./materialProperty Av = getattr(self, 'ave'+AvType+'2CC') - V = Utils.sdiag(self.vol) - return self.dim * Utils.sdiag(Av.T * V * materialProperty) + Vprop = self.vol * Utils.mkvc(materialProperty) + return self.dim * Utils.sdiag(Av.T * Vprop) if materialProperty.size == self.nC*self.dim: - if invertProperty: - materialProperty = 1./materialProperty Av = getattr(self, 'ave'+AvType+'2CCV') V = sp.kron(sp.identity(self.dim), Utils.sdiag(self.vol)) return Utils.sdiag(Av.T * V * Utils.mkvc(materialProperty)) @@ -576,11 +576,20 @@ class TensorMesh(BaseRectangularMesh, TensorView, DiffOperators, InnerProducts): :rtype: scipy.csr_matrix :return: M, the inner product matrix (nF, nF) """ - if materialProperty is None or materialProperty.size == self.nC: + if materialProperty is None: + return None + if Utils.isScalar(materialProperty): + Av = getattr(self, 'ave'+AvType+'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 materialProperty.size == self.nC: Av = getattr(self, 'ave'+AvType+'2CC') V = Utils.sdiag(self.vol) if v is None: - return self.dim * Av.T * Utils.sdiag(self.vol) + return self.dim * Av.T * V return Utils.sdiag(v) * self.dim * Av.T * V if materialProperty.size == self.nC*self.dim: # anisotropic Av = getattr(self, 'ave'+AvType+'2CCV') diff --git a/SimPEG/Tests/test_massMatrices.py b/SimPEG/Tests/test_massMatrices.py index 3d3946ad..67de753e 100644 --- a/SimPEG/Tests/test_massMatrices.py +++ b/SimPEG/Tests/test_massMatrices.py @@ -122,7 +122,7 @@ class TestInnerProducts2D(OrderTest): sigma = np.c_[call(sigma1, Gc), call(sigma2, Gc)] analytic = 189959./120 # Found using sympy. z=5 elif self.sigmaTest == 3: - sigma = np.c_[call(sigma1, Gc), call(sigma2, Gc), call(sigma3, Gc)] + sigma = np.r_[call(sigma1, Gc), call(sigma2, Gc), call(sigma3, Gc)] analytic = 781427./360 # Found using sympy. z=5 if self.location == 'edges': diff --git a/SimPEG/Tests/test_massMatricesDerivs.py b/SimPEG/Tests/test_massMatricesDerivs.py index 8be5e4eb..d66c0b65 100644 --- a/SimPEG/Tests/test_massMatricesDerivs.py +++ b/SimPEG/Tests/test_massMatricesDerivs.py @@ -16,7 +16,7 @@ class TestInnerProductsDerivs(unittest.TestCase): return M*v, Md Md = mesh.getFaceInnerProductDeriv(sig, doFast=fast) return M*v, Utils.sdiag(v)*Md - sig = np.random.rand(mesh.nC*rep) + sig = np.random.rand(1) if rep is 0 else np.random.rand(mesh.nC*rep) return checkDerivative(fun, sig, num=5, plotIt=False) def doTestEdge(self, h, rep, vec, fast): @@ -29,9 +29,15 @@ class TestInnerProductsDerivs(unittest.TestCase): return M*v, Md Md = mesh.getEdgeInnerProductDeriv(sig, doFast=fast) return M*v, Utils.sdiag(v)*Md - sig = np.random.rand(mesh.nC*rep) + sig = np.random.rand(1) if rep is 0 else np.random.rand(mesh.nC*rep) return checkDerivative(fun, sig, num=5, plotIt=False) + def test_FaceIP_1D_float(self): + self.assertTrue(self.doTestFace([10],0,True, False)) + def test_FaceIP_2D_float(self): + self.assertTrue(self.doTestFace([10, 4],0,True, False)) + def test_FaceIP_3D_float(self): + self.assertTrue(self.doTestFace([10, 4, 5],0,True, False)) def test_FaceIP_1D_isotropic(self): self.assertTrue(self.doTestFace([10],1,True, False)) def test_FaceIP_2D_isotropic(self): @@ -47,6 +53,12 @@ class TestInnerProductsDerivs(unittest.TestCase): def test_FaceIP_3D_tensor(self): self.assertTrue(self.doTestFace([10, 4, 5],6,True, False)) + def test_FaceIP_1D_float_fast(self): + self.assertTrue(self.doTestFace([10],0, False, True)) + def test_FaceIP_2D_float_fast(self): + self.assertTrue(self.doTestFace([10, 4],0, False, True)) + def test_FaceIP_3D_float_fast(self): + self.assertTrue(self.doTestFace([10, 4, 5],0, False, True)) def test_FaceIP_1D_isotropic_fast(self): self.assertTrue(self.doTestFace([10],1, False, True)) def test_FaceIP_2D_isotropic_fast(self): @@ -59,6 +71,10 @@ class TestInnerProductsDerivs(unittest.TestCase): self.assertTrue(self.doTestFace([10, 4, 5],3, False, True)) + def test_EdgeIP_2D_float(self): + self.assertTrue(self.doTestEdge([10, 4],0,True, False)) + def test_EdgeIP_3D_float(self): + self.assertTrue(self.doTestEdge([10, 4, 5],0,True, False)) def test_EdgeIP_2D_isotropic(self): self.assertTrue(self.doTestEdge([10, 4],1,True, False)) def test_EdgeIP_3D_isotropic(self): @@ -72,6 +88,10 @@ class TestInnerProductsDerivs(unittest.TestCase): def test_EdgeIP_3D_tensor(self): self.assertTrue(self.doTestEdge([10, 4, 5],6,True, False)) + def test_EdgeIP_2D_float_fast(self): + self.assertTrue(self.doTestEdge([10, 4],0, False, True)) + def test_EdgeIP_3D_float_fast(self): + self.assertTrue(self.doTestEdge([10, 4, 5],0, False, True)) def test_EdgeIP_2D_isotropic_fast(self): self.assertTrue(self.doTestEdge([10, 4],1, False, True)) def test_EdgeIP_3D_isotropic_fast(self): diff --git a/SimPEG/Tests/test_utils.py b/SimPEG/Tests/test_utils.py index e69c8473..ad737489 100644 --- a/SimPEG/Tests/test_utils.py +++ b/SimPEG/Tests/test_utils.py @@ -164,6 +164,12 @@ class TestSequenceFunctions(unittest.TestCase): Z = B2*A - sp.identity(M.nC*3) self.assertTrue(np.linalg.norm(Z.todense().ravel(), 2) < TOL) + def test_isFloat(self): + self.assertTrue(isScalar(1.)) + self.assertTrue(isScalar(1)) + self.assertTrue(isScalar(long(1))) + self.assertTrue(isScalar(np.r_[1.])) + self.assertTrue(isScalar(np.r_[1])) if __name__ == '__main__': unittest.main() diff --git a/SimPEG/Utils/matutils.py b/SimPEG/Utils/matutils.py index c15f6234..f1ea7ec7 100644 --- a/SimPEG/Utils/matutils.py +++ b/SimPEG/Utils/matutils.py @@ -1,6 +1,15 @@ import numpy as np import scipy.sparse as sp + +def isScalar(f): + scalarTypes = [float, int, long, np.float_, np.int_] + if type(f) in scalarTypes: + return True + elif type(f) == np.ndarray and f.size == 1 and type(f[0]) in scalarTypes: + return True + return False + def mkvc(x, numDims=1): """Creates a vector with the number of dimension specified @@ -248,7 +257,7 @@ def makePropertyTensor(M, sigma): if sigma is None: # default is ones sigma = np.ones(M.nC) - if type(sigma) in [float, int, long]: + if isScalar(sigma): sigma = sigma * np.ones(M.nC) if M.dim == 1: @@ -294,7 +303,7 @@ def invPropertyTensor(M, tensor, returnMatrix=False): T = None - if type(tensor) in [float, int, long]: + if isScalar(tensor): T = 1./tensor elif tensor.size == M.nC: # Isotropic!