diff --git a/SimPEG/Mesh/InnerProducts.py b/SimPEG/Mesh/InnerProducts.py index 9bc092d3..f0e45665 100644 --- a/SimPEG/Mesh/InnerProducts.py +++ b/SimPEG/Mesh/InnerProducts.py @@ -148,18 +148,16 @@ class InnerProducts(object): :rtype: scipy.csr_matrix :return: dMdm, the derivative of the inner product matrix (nE, nC*nA) """ - tensorType = TensorType(self, prop) fast = None - if hasattr(self, '_fastInnerProductDeriv') and doFast: - fast = self._fastInnerProductDeriv(projType, tensorType, invProp=invProp, invMat=invMat) + fast = self._fastInnerProductDeriv(projType, prop, invProp=invProp, invMat=invMat) + if fast is not None: + return fast if invProp or invMat: raise NotImplementedError('inverting the property or the matrix is not yet implemented for this mesh/tensorType. You should write it!') - if fast is not None: - return fast - + tensorType = TensorType(self, prop) P = self._getInnerProductProjectionMatrices(projType, tensorType=tensorType) def innerProductDeriv(v): return self._getInnerProductDerivFunction(tensorType, P, projType, v) diff --git a/SimPEG/Mesh/TensorMesh.py b/SimPEG/Mesh/TensorMesh.py index 854a111a..d9a14aaa 100644 --- a/SimPEG/Mesh/TensorMesh.py +++ b/SimPEG/Mesh/TensorMesh.py @@ -300,46 +300,57 @@ class BaseTensorMesh(BaseRectangularMesh): else: return M - def _fastInnerProductDeriv(self, projType, tensorType): + def _fastInnerProductDeriv(self, projType, prop, invProp=False, invMat=False): """ :param str projType: 'E' or 'F' :param TensorType tensorType: type of the tensor + :param bool invProp: inverts the material property + :param bool invMat: inverts the matrix :rtype: function :return: dMdmu, the derivative of the inner product matrix """ assert projType in ['F', 'E'], "projType must be 'F' for faces or 'E' for edges" - if tensorType == -1: - return None + tensorType = Utils.TensorType(self, prop) + + dMdprop = None + + if invMat: + MI = self._fastInnerProduct(projType, prop, invProp=invProp, invMat=invMat) 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 - dim_x_AvT_x_V_x_ones = self.dim * Av.T * V * ones - def scalarInnerProductDeriv(v): - return Utils.sdiag(v) * dim_x_AvT_x_V_x_ones - return scalarInnerProductDeriv + if not invMat and not invProp: + dMdprop = self.dim * Av.T * V * ones + elif invMat and invProp: + dMdprop = self.dim * Utils.sdiag(MI.diagonal()**2) * Av.T * V * ones * Utils.sdiag(1./prop**2) if tensorType == 1: Av = getattr(self, 'ave'+projType+'2CC') V = Utils.sdiag(self.vol) - dim_x_AvT_x_V = self.dim * Av.T * V - def isotropicInnerProductDeriv(v=None): - if v is None: - print 'Depreciation Warning: TensorMesh.isotropicInnerProductDeriv. You should be supplying a vector. Returning: d * Av.T * V' - return dim_x_AvT_x_V - return Utils.sdiag(v) * dim_x_AvT_x_V - return isotropicInnerProductDeriv + if not invMat and not invProp: + dMdprop = self.dim * Av.T * V + elif invMat and invProp: + dMdprop = self.dim * Utils.sdiag(MI.diagonal()**2) * Av.T * V * Utils.sdiag(1./prop**2) if tensorType == 2: # anisotropic Av = getattr(self, 'ave'+projType+'2CCV') V = sp.kron(sp.identity(self.dim), Utils.sdiag(self.vol)) - AvT_x_V = Av.T * V - def anisotropicInnerProductDeriv(v): - return Utils.sdiag(v) * AvT_x_V - return anisotropicInnerProductDeriv + if not invMat and not invProp: + dMdprop = Av.T * V + elif invMat and invProp: + dMdprop = Utils.sdiag(MI.diagonal()**2) * Av.T * V * Utils.sdiag(1./prop**2) + + if dMdprop is not None: + def innerProductDeriv(v=None): + if v is None: + print 'Depreciation Warning: TensorMesh.innerProductDeriv. You should be supplying a vector. Use: sdiag(u)*dMdprop' + return dMdprop + return Utils.sdiag(v) * dMdprop + return innerProductDeriv + else: + return None diff --git a/SimPEG/Tests/test_massMatricesDerivs.py b/SimPEG/Tests/test_massMatricesDerivs.py index 5eac4ca8..4b7d63cd 100644 --- a/SimPEG/Tests/test_massMatricesDerivs.py +++ b/SimPEG/Tests/test_massMatricesDerivs.py @@ -20,7 +20,7 @@ class TestInnerProductsDerivs(unittest.TestCase): M = mesh.getFaceInnerProduct(sig, invProp=invProp, invMat=invMat) Md = mesh.getFaceInnerProductDeriv(sig, invProp=invProp, invMat=invMat, doFast=fast) return M*v, Md(v) - print meshType, 'Face', h, rep, fast + print meshType, 'Face', h, rep, fast, ('harmonic' if invProp and invMat else 'standard') return checkDerivative(fun, sig, num=5, plotIt=False) def doTestEdge(self, h, rep, fast, meshType, invProp=False, invMat=False): @@ -37,7 +37,7 @@ class TestInnerProductsDerivs(unittest.TestCase): M = mesh.getEdgeInnerProduct(sig, invProp=invProp, invMat=invMat) Md = mesh.getEdgeInnerProductDeriv(sig, invProp=invProp, invMat=invMat, doFast=fast) return M*v, Md(v) - print meshType, 'Edge', h, rep, fast + print meshType, 'Edge', h, rep, fast, ('harmonic' if invProp and invMat else 'standard') return checkDerivative(fun, sig, num=5, plotIt=False) def test_FaceIP_1D_float(self): @@ -118,22 +118,22 @@ class TestInnerProductsDerivs(unittest.TestCase): - # def test_FaceIP_1D_float_fast_harmonic(self): - # self.assertTrue(self.doTestFace([10],0, True, 'Tensor', invProp=True, invMat=True)) - # def test_FaceIP_2D_float_fast_harmonic(self): - # self.assertTrue(self.doTestFace([10, 4],0, True, 'Tensor', invProp=True, invMat=True)) - # def test_FaceIP_3D_float_fast_harmonic(self): - # self.assertTrue(self.doTestFace([10, 4, 5],0, True, 'Tensor', invProp=True, invMat=True)) + def test_FaceIP_1D_float_fast_harmonic(self): + self.assertTrue(self.doTestFace([10],0, True, 'Tensor', invProp=True, invMat=True)) + def test_FaceIP_2D_float_fast_harmonic(self): + self.assertTrue(self.doTestFace([10, 4],0, True, 'Tensor', invProp=True, invMat=True)) + def test_FaceIP_3D_float_fast_harmonic(self): + self.assertTrue(self.doTestFace([10, 4, 5],0, True, 'Tensor', invProp=True, invMat=True)) def test_FaceIP_1D_isotropic_fast_harmonic(self): self.assertTrue(self.doTestFace([10],1, True, 'Tensor', invProp=True, invMat=True)) def test_FaceIP_2D_isotropic_fast_harmonic(self): self.assertTrue(self.doTestFace([10, 4],1, True, 'Tensor', invProp=True, invMat=True)) def test_FaceIP_3D_isotropic_fast_harmonic(self): self.assertTrue(self.doTestFace([10, 4, 5],1, True, 'Tensor', invProp=True, invMat=True)) - # def test_FaceIP_2D_anisotropic_fast_harmonic(self): - # self.assertTrue(self.doTestFace([10, 4],2, True, 'Tensor', invProp=True, invMat=True)) - # def test_FaceIP_3D_anisotropic_fast_harmonic(self): - # self.assertTrue(self.doTestFace([10, 4, 5],3, True, 'Tensor', invProp=True, invMat=True)) + def test_FaceIP_2D_anisotropic_fast_harmonic(self): + self.assertTrue(self.doTestFace([10, 4],2, True, 'Tensor', invProp=True, invMat=True)) + def test_FaceIP_3D_anisotropic_fast_harmonic(self): + self.assertTrue(self.doTestFace([10, 4, 5],3, True, 'Tensor', invProp=True, invMat=True))