diff --git a/SimPEG/Mesh/InnerProducts.py b/SimPEG/Mesh/InnerProducts.py index 1adbb796..8268b29a 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, isScalar +from SimPEG.Utils import * import numpy as np @@ -11,11 +11,12 @@ class InnerProducts(object): 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, doFast=True): + 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) @@ -23,7 +24,7 @@ class InnerProducts(object): fast = None if returnP is False and hasattr(self, '_fastFaceInnerProduct') and doFast: - fast = self._fastFaceInnerProduct(prop=prop, invProp=invProp) + fast = self._fastFaceInnerProduct(prop=prop, invProp=invProp, invMat=invMat) if fast is not None: return fast @@ -67,6 +68,12 @@ class InnerProducts(object): 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: @@ -95,11 +102,12 @@ class InnerProducts(object): return self._getInnerProductDeriv(prop, v, P, self.nF) def getEdgeInnerProduct(self, prop=None, returnP=False, - invProp=False, doFast=True): + 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) @@ -107,7 +115,7 @@ class InnerProducts(object): fast = None if returnP is False and hasattr(self, '_fastEdgeInnerProduct') and doFast: - fast = self._fastEdgeInnerProduct(prop=prop, invProp=invProp) + fast = self._fastEdgeInnerProduct(prop=prop, invProp=invProp, invMat=invMat) if fast is not None: return fast @@ -146,6 +154,12 @@ class InnerProducts(object): 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 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: diff --git a/SimPEG/Mesh/TensorMesh.py b/SimPEG/Mesh/TensorMesh.py index 2f9a65a6..72166b1f 100644 --- a/SimPEG/Mesh/TensorMesh.py +++ b/SimPEG/Mesh/TensorMesh.py @@ -241,7 +241,7 @@ class BaseTensorMesh(BaseRectangularMesh): return Q.tocsr() - def _fastFaceInnerProduct(self, prop=None, invProp=False): + def _fastFaceInnerProduct(self, prop=None, invProp=False, invMat=False): """ Fast version of getFaceInnerProduct. This does not handle the case of a full tensor prop. @@ -249,13 +249,14 @@ class BaseTensorMesh(BaseRectangularMesh): :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 (nF, nF) """ - return self._fastInnerProduct('F', prop=prop, invProp=invProp) + return self._fastInnerProduct('F', prop=prop, invProp=invProp, invMat=invMat) - def _fastEdgeInnerProduct(self, prop=None, invProp=False): + def _fastEdgeInnerProduct(self, prop=None, invProp=False, invMat=False): """ Fast version of getEdgeInnerProduct. This does not handle the case of a full tensor prop. @@ -263,13 +264,14 @@ class BaseTensorMesh(BaseRectangularMesh): :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) + return self._fastInnerProduct('E', prop=prop, invProp=invProp, invMat=invMat) - def _fastInnerProduct(self, AvType, prop=None, invProp=False): + 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. @@ -278,6 +280,7 @@ class BaseTensorMesh(BaseRectangularMesh): :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) """ @@ -293,12 +296,18 @@ class BaseTensorMesh(BaseRectangularMesh): if prop.size == self.nC: Av = getattr(self, 'ave'+AvType+'2CC') Vprop = self.vol * Utils.mkvc(prop) - return self.dim * Utils.sdiag(Av.T * Vprop) - if prop.size == self.nC*self.dim: + M = self.dim * Utils.sdiag(Av.T * Vprop) + elif prop.size == self.nC*self.dim: 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(prop)) + M = Utils.sdiag(Av.T * V * Utils.mkvc(prop)) + else: + return None + if invMat: + return Utils.sdInv(M) + else: + return M def _fastFaceInnerProductDeriv(self, prop=None, v=None): """ diff --git a/SimPEG/Tests/test_massMatrices.py b/SimPEG/Tests/test_massMatrices.py index 67de753e..b167ce9e 100644 --- a/SimPEG/Tests/test_massMatrices.py +++ b/SimPEG/Tests/test_massMatrices.py @@ -1,6 +1,7 @@ import numpy as np import unittest from TestUtils import OrderTest +from SimPEG import Utils class TestInnerProducts(OrderTest): @@ -130,14 +131,21 @@ class TestInnerProducts2D(OrderTest): Ec = np.vstack((cart(self.M.gridEx), cart(self.M.gridEy))) E = self.M.projectEdgeVector(Ec) - A = self.M.getEdgeInnerProduct(sigma) + if self.invProp: + A = self.M.getEdgeInnerProduct(Utils.invPropertyTensor(self.M, sigma), invProp=True) + else: + A = self.M.getEdgeInnerProduct(sigma) numeric = E.T.dot(A.dot(E)) elif self.location == 'faces': cart = lambda g: np.c_[call(ex, g), call(ey, g)] Fc = np.vstack((cart(self.M.gridFx), cart(self.M.gridFy))) F = self.M.projectFaceVector(Fc) - A = self.M.getFaceInnerProduct(sigma) + + if self.invProp: + A = self.M.getFaceInnerProduct(Utils.invPropertyTensor(self.M, sigma), invProp=True) + else: + A = self.M.getFaceInnerProduct(sigma) numeric = F.T.dot(A.dot(F)) err = np.abs(numeric - analytic) @@ -147,36 +155,60 @@ class TestInnerProducts2D(OrderTest): self.name = "2D Edge Inner Product - Isotropic" self.location = 'edges' self.sigmaTest = 1 + self.invProp = True + self.orderTest() + self.name += " - invProp" + self.invProp = False self.orderTest() def test_order3_edges(self): self.name = "2D Edge Inner Product - Anisotropic" self.location = 'edges' self.sigmaTest = 2 + self.invProp = True + self.orderTest() + self.name += " - invProp" + self.invProp = False self.orderTest() def test_order6_edges(self): self.name = "2D Edge Inner Product - Full Tensor" self.location = 'edges' self.sigmaTest = 3 + self.invProp = True + self.orderTest() + self.name += " - invProp" + self.invProp = False self.orderTest() def test_order1_faces(self): self.name = "2D Face Inner Product - Isotropic" self.location = 'faces' self.sigmaTest = 1 + self.invProp = True + self.orderTest() + self.name += " - invProp" + self.invProp = False self.orderTest() def test_order2_faces(self): self.name = "2D Face Inner Product - Anisotropic" self.location = 'faces' self.sigmaTest = 2 + self.invProp = True + self.orderTest() + self.name += " - invProp" + self.invProp = False self.orderTest() def test_order3_faces(self): self.name = "2D Face Inner Product - Full Tensor" self.location = 'faces' self.sigmaTest = 3 + self.invProp = True + self.orderTest() + self.name += " - invProp" + self.invProp = False self.orderTest() @@ -205,7 +237,10 @@ class TestInnerProducts1D(OrderTest): if self.location == 'faces': F = call(ex, self.M.gridFx) - A = self.M.getFaceInnerProduct(sigma) + if self.invProp: + A = self.M.getFaceInnerProduct(1/sigma, invProp=True) + else: + A = self.M.getFaceInnerProduct(sigma) numeric = F.T.dot(A.dot(F)) err = np.abs(numeric - analytic) @@ -215,6 +250,10 @@ class TestInnerProducts1D(OrderTest): self.name = "1D Face Inner Product" self.location = 'faces' self.sigmaTest = 1 + self.invProp = True + self.orderTest() + self.name += " - invProp" + self.invProp = False self.orderTest() diff --git a/SimPEG/Tests/test_massMatricesDerivs.py b/SimPEG/Tests/test_massMatricesDerivs.py index d66c0b65..1da88173 100644 --- a/SimPEG/Tests/test_massMatricesDerivs.py +++ b/SimPEG/Tests/test_massMatricesDerivs.py @@ -70,6 +70,22 @@ class TestInnerProductsDerivs(unittest.TestCase): 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)) def test_EdgeIP_2D_float(self): self.assertTrue(self.doTestEdge([10, 4],0,True, False)) @@ -101,6 +117,18 @@ class TestInnerProductsDerivs(unittest.TestCase): 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)) diff --git a/SimPEG/Tests/test_utils.py b/SimPEG/Tests/test_utils.py index 0ac03ef3..4ffb66a3 100644 --- a/SimPEG/Tests/test_utils.py +++ b/SimPEG/Tests/test_utils.py @@ -140,6 +140,39 @@ 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): + M = Mesh.TensorMesh([6, 6]) + a1 = np.random.rand(M.nC) + a2 = np.random.rand(M.nC) + a3 = np.random.rand(M.nC) + prop1 = a1 + prop2 = np.c_[a1, a2] + prop3 = np.c_[a1, a2, a3] + + for ii, prop in enumerate([4, prop1, prop2, prop3]): + self.assertTrue(tensorType(M, prop) == ii) + + self.assertRaises(Exception, tensorType, M, np.c_[a1, a2, a3, a3]) + self.assertTrue(tensorType(M, None) == -1) + + def test_tensorType3D(self): + M = Mesh.TensorMesh([6, 6, 7]) + a1 = np.random.rand(M.nC) + a2 = np.random.rand(M.nC) + a3 = np.random.rand(M.nC) + a4 = np.random.rand(M.nC) + a5 = np.random.rand(M.nC) + a6 = np.random.rand(M.nC) + prop1 = a1 + prop2 = np.c_[a1, a2, a3] + 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.assertRaises(Exception, tensorType, M, np.c_[a1, a2, a3, a3]) + self.assertTrue(tensorType(M, None) == -1) + def test_invPropertyTensor3D(self): M = Mesh.TensorMesh([6, 6, 6]) diff --git a/SimPEG/Utils/matutils.py b/SimPEG/Utils/matutils.py index 992ed3bc..2d3b87f7 100644 --- a/SimPEG/Utils/matutils.py +++ b/SimPEG/Utils/matutils.py @@ -251,85 +251,79 @@ def inv2X2BlockDiagonal(a11, a12, a21, a22, returnMatrix=True): return sp.vstack((sp.hstack((sdiag(b11), sdiag(b12))), sp.hstack((sdiag(b21), sdiag(b22))))) -def makePropertyTensor(M, sigma): - if sigma is None: # default is ones - sigma = np.ones(M.nC) +def tensorType(M, tensor): + if tensor is None: # default is ones + return -1 - if isScalar(sigma): - sigma = sigma * np.ones(M.nC) + if isScalar(tensor): + return 0 - if M.dim == 1: - if sigma.size == M.nC: # Isotropic! - sigma = mkvc(sigma) # ensure it is a vector. - Sigma = sdiag(sigma) - else: - raise Exception('Unexpected shape of sigma') - elif M.dim == 2: - if sigma.size == M.nC: # Isotropic! - sigma = mkvc(sigma) # ensure it is a vector. - Sigma = sdiag(np.r_[sigma, sigma]) - elif sigma.size == M.nC*2: # Diagonal tensor - sigma = sigma.reshape((M.nC,2), order='F') - Sigma = sdiag(np.r_[sigma[:, 0], sigma[:, 1]]) - elif sigma.size == M.nC*3: # Fully anisotropic - sigma = sigma.reshape((M.nC,3), order='F') - row1 = sp.hstack((sdiag(sigma[:, 0]), sdiag(sigma[:, 2]))) - row2 = sp.hstack((sdiag(sigma[:, 2]), sdiag(sigma[:, 1]))) - Sigma = sp.vstack((row1, row2)) - else: - raise Exception('Unexpected shape of sigma') - elif M.dim == 3: - if sigma.size == M.nC: # Isotropic! - sigma = mkvc(sigma) # ensure it is a vector. - Sigma = sdiag(np.r_[sigma, sigma, sigma]) - elif sigma.size == M.nC*3: # Diagonal tensor - sigma = sigma.reshape((M.nC,3), order='F') - Sigma = sdiag(np.r_[sigma[:, 0], sigma[:, 1], sigma[:, 2]]) - elif sigma.size == M.nC*6: # Fully anisotropic - sigma = sigma.reshape((M.nC,6), order='F') - row1 = sp.hstack((sdiag(sigma[:, 0]), sdiag(sigma[:, 3]), sdiag(sigma[:, 4]))) - row2 = sp.hstack((sdiag(sigma[:, 3]), sdiag(sigma[:, 1]), sdiag(sigma[:, 5]))) - row3 = sp.hstack((sdiag(sigma[:, 4]), sdiag(sigma[:, 5]), sdiag(sigma[:, 2]))) - Sigma = sp.vstack((row1, row2, row3)) - else: - raise Exception('Unexpected shape of sigma') + 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') + +def makePropertyTensor(M, tensor): + if tensor is None: # default is ones + tensor = np.ones(M.nC) + + if isScalar(tensor): + tensor = tensor * np.ones(M.nC) + + propType = tensorType(M, tensor) + if propType == 1: # Isotropic! + Sigma = sp.kron(sp.identity(M.dim), sdiag(mkvc(tensor))) + elif propType == 2: # Diagonal tensor + Sigma = sdiag(mkvc(tensor)) + elif M.dim == 2 and tensor.size == M.nC*3: # Fully anisotropic, 2D + tensor = tensor.reshape((M.nC,3), order='F') + row1 = sp.hstack((sdiag(tensor[:, 0]), sdiag(tensor[:, 2]))) + row2 = sp.hstack((sdiag(tensor[:, 2]), sdiag(tensor[:, 1]))) + Sigma = sp.vstack((row1, row2)) + elif M.dim == 3 and tensor.size == M.nC*6: # Fully anisotropic, 3D + tensor = tensor.reshape((M.nC,6), order='F') + row1 = sp.hstack((sdiag(tensor[:, 0]), sdiag(tensor[:, 3]), sdiag(tensor[:, 4]))) + row2 = sp.hstack((sdiag(tensor[:, 3]), sdiag(tensor[:, 1]), sdiag(tensor[:, 5]))) + row3 = sp.hstack((sdiag(tensor[:, 4]), sdiag(tensor[:, 5]), sdiag(tensor[:, 2]))) + Sigma = sp.vstack((row1, row2, row3)) + else: + raise Exception('Unexpected shape of tensor') return Sigma def invPropertyTensor(M, tensor, returnMatrix=False): - T = None + propType = tensorType(M, tensor) if isScalar(tensor): T = 1./tensor - - elif tensor.size == M.nC: # Isotropic! + elif propType < 3: # Isotropic or Diagonal T = 1./mkvc(tensor) # ensure it is a vector. - - elif M.dim == 2: - if tensor.size == M.nC*2: # Diagonal tensor - T = 1./tensor - elif tensor.size == M.nC*3: # Fully anisotropic - tensor = tensor.reshape((M.nC,3), order='F') - B = inv2X2BlockDiagonal(tensor[:,0], tensor[:,2], - tensor[:,2], tensor[:,1], - returnMatrix=False) - b11, b12, b21, b22 = B - T = np.r_[b11, b22, b12] - elif M.dim == 3: - if tensor.size == M.nC*3: # Diagonal tensor - T = 1./tensor - elif tensor.size == M.nC*6: # Fully anisotropic - tensor = tensor.reshape((M.nC,6), order='F') - B = inv3X3BlockDiagonal(tensor[:,0], tensor[:,3], tensor[:,4], - tensor[:,3], tensor[:,1], tensor[:,5], - tensor[:,4], tensor[:,5], tensor[:,2], - returnMatrix=False) - b11, b12, b13, b21, b22, b23, b31, b32, b33 = B - T = np.r_[b11, b22, b33, b12, b13, b23] - - if T is None: + elif M.dim == 2 and tensor.size == M.nC*3: # Fully anisotropic, 2D + tensor = tensor.reshape((M.nC,3), order='F') + B = inv2X2BlockDiagonal(tensor[:,0], tensor[:,2], + tensor[:,2], tensor[:,1], + returnMatrix=False) + b11, b12, b21, b22 = B + T = np.r_[b11, b22, b12] + elif M.dim == 3 and tensor.size == M.nC*6: # Fully anisotropic, 3D + tensor = tensor.reshape((M.nC,6), order='F') + B = inv3X3BlockDiagonal(tensor[:,0], tensor[:,3], tensor[:,4], + tensor[:,3], tensor[:,1], tensor[:,5], + tensor[:,4], tensor[:,5], tensor[:,2], + returnMatrix=False) + b11, b12, b13, b21, b22, b23, b31, b32, b33 = B + T = np.r_[b11, b22, b33, b12, b13, b23] + else: raise Exception('Unexpected shape of tensor') if returnMatrix: