From 3dd54bc0e3733301e69d54d89c9296fafcf98902 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Mon, 3 Mar 2014 11:01:12 -0800 Subject: [PATCH] Full anisotropic tensor derivatives. --- SimPEG/Mesh/InnerProducts.py | 123 ++++++++++++++++++++---- SimPEG/Mesh/TensorMesh.py | 19 ++-- SimPEG/Tests/test_massMatricesDerivs.py | 107 ++++++++++++++------- 3 files changed, 186 insertions(+), 63 deletions(-) diff --git a/SimPEG/Mesh/InnerProducts.py b/SimPEG/Mesh/InnerProducts.py index 266661a7..bb77d5e5 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 +from SimPEG.Utils import sub2ind, ndgrid, mkvc, getSubArray, sdiag, inv3X3BlockDiagonal, inv2X2BlockDiagonal, makePropertyTensor, invPropertyTensor, spzeros import numpy as np @@ -11,7 +11,7 @@ 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, materialProperty=None, returnP=False, - invertProperty=False): + invertProperty=False, doFast=True): """ :param numpy.array materialProperty: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6)) :param bool returnP: returns the projection matrices @@ -21,7 +21,7 @@ class InnerProducts(object): """ fast = None - if returnP is False and hasattr(self, '_fastFaceInnerProduct'): + if returnP is False and hasattr(self, '_fastFaceInnerProduct') and doFast: fast = self._fastFaceInnerProduct(materialProperty=materialProperty, invertProperty=invertProperty) if fast is not None: @@ -71,7 +71,7 @@ class InnerProducts(object): else: return A - def getFaceInnerProductDeriv(self, materialProperty=None, P=None): + def getFaceInnerProductDeriv(self, materialProperty=None, v=None, P=None, doFast=True): """ :param numpy.array materialProperty: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6)) :rtype: scipy.csr_matrix @@ -79,24 +79,18 @@ class InnerProducts(object): """ fast = None - if hasattr(self, '_fastFaceInnerProductDeriv'): - fast = self._fastFaceInnerProductDeriv(materialProperty=materialProperty) + if hasattr(self, '_fastFaceInnerProductDeriv') and doFast: + fast = self._fastFaceInnerProductDeriv(materialProperty=materialProperty, v=v) if fast is not None: return fast - raise NotImplementedError('Derivatives for the material property specified are not yet implemented.') - if P is None: - M, P = getFaceInnerProduct(self, materialProperty=materialProperty, returnP=True) + M, P = self.getFaceInnerProduct(materialProperty=materialProperty, returnP=True) - d = self.dim + return self._getInnerProductDeriv(materialProperty, v, P, self.nF) - if d == 1: - P[0].T * sp.identity(n) * P[0] - - - def getEdgeInnerProduct(self, materialProperty=None, returnP=False, invertProperty=False): + def getEdgeInnerProduct(self, materialProperty=None, returnP=False, invertProperty=False, doFast=True): """ :param numpy.array materialProperty: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6)) :param bool returnP: returns the projection matrices @@ -106,7 +100,7 @@ class InnerProducts(object): """ fast = None - if returnP is False and hasattr(self, '_fastEdgeInnerProduct'): + if returnP is False and hasattr(self, '_fastEdgeInnerProduct') and doFast: fast = self._fastEdgeInnerProduct(materialProperty=materialProperty, invertProperty=invertProperty) if fast is not None: @@ -151,8 +145,7 @@ class InnerProducts(object): else: return A - - def getEdgeInnerProductDeriv(self, materialProperty=None, P=None): + def getEdgeInnerProductDeriv(self, materialProperty=None, v=None, P=None, doFast=True): """ :param numpy.array materialProperty: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6)) :rtype: scipy.csr_matrix @@ -161,14 +154,102 @@ class InnerProducts(object): fast = None - if hasattr(self, '_fastEdgeInnerProductDeriv'): - fast = self._fastEdgeInnerProductDeriv(materialProperty=materialProperty) + if hasattr(self, '_fastEdgeInnerProductDeriv') and doFast: + fast = self._fastEdgeInnerProductDeriv(materialProperty=materialProperty, v=v) if fast is not None: return fast - raise NotImplementedError('Derivatives for the material property specified are not yet implemented.') + if P is None: + M, P = self.getEdgeInnerProduct(materialProperty=materialProperty, returnP=True) + return self._getInnerProductDeriv(materialProperty, v, P, self.nE) + + def _getInnerProductDeriv(self, materialProperty, v, P, n): + + 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) + 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: + 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: + dMdm1 = spzeros(n, self.nC) + dMdm2 = spzeros(n, self.nC) + for i, p in enumerate(P): + Y = p * v + y1 = Y[:self.nC] + y2 = Y[self.nC:] + dMdm1 = dMdm1 + p.T * sp.vstack(( sdiag( y1 ), Z)) + dMdm2 = dMdm2 + p.T * sp.vstack(( Z, sdiag( y2 ))) + dMdm = sp.hstack((dMdm1, dMdm2)) + if materialProperty.size == self.nC*3: + dMdm1 = spzeros(n, self.nC) + dMdm2 = spzeros(n, self.nC) + dMdm3 = spzeros(n, self.nC) + for i, p in enumerate(P): + Y = p * v + y1 = Y[:self.nC] + y2 = Y[self.nC:] + dMdm1 = dMdm1 + p.T * sp.vstack(( sdiag( y1 ), Z)) + dMdm2 = dMdm2 + p.T * sp.vstack(( Z, sdiag( y2 ))) + dMdm3 = dMdm3 + p.T * sp.vstack(( sdiag( y2 ), sdiag( y1 ))) + dMdm = sp.hstack((dMdm1, dMdm2, dMdm3)) + elif d == 3: + if materialProperty is None or 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: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: + dMdm1 = spzeros(n, self.nC) + dMdm2 = spzeros(n, self.nC) + dMdm3 = spzeros(n, self.nC) + for i, p in enumerate(P): + Y = p * v + y1 = Y[:self.nC] + y2 = Y[self.nC:self.nC*2] + y3 = Y[self.nC*2:] + dMdm1 = dMdm1 + p.T * sp.vstack(( sdiag( y1 ), Z, Z)) + dMdm2 = dMdm2 + p.T * sp.vstack(( Z, sdiag( y2 ), Z)) + dMdm3 = dMdm3 + p.T * sp.vstack(( Z, Z, sdiag( y3 ))) + dMdm = sp.hstack((dMdm1, dMdm2, dMdm3)) + if materialProperty.size == self.nC*6: + dMdm1 = spzeros(n, self.nC) + dMdm2 = spzeros(n, self.nC) + dMdm3 = spzeros(n, self.nC) + dMdm4 = spzeros(n, self.nC) + dMdm5 = spzeros(n, self.nC) + dMdm6 = spzeros(n, self.nC) + for i, p in enumerate(P): + Y = p * v + y1 = Y[:self.nC] + y2 = Y[self.nC:self.nC*2] + y3 = Y[self.nC*2:] + dMdm1 = dMdm1 + p.T * sp.vstack(( sdiag( y1 ), Z, Z)) + dMdm2 = dMdm2 + p.T * sp.vstack(( Z, sdiag( y2 ), Z)) + dMdm3 = dMdm3 + p.T * sp.vstack(( Z, Z, sdiag( y3 ))) + dMdm4 = dMdm4 + p.T * sp.vstack(( sdiag( y2 ), sdiag( y1 ), Z)) + dMdm5 = dMdm5 + p.T * sp.vstack(( sdiag( y3 ), Z, sdiag( y1 ))) + dMdm6 = dMdm6 + p.T * sp.vstack(( Z, sdiag( y3 ), sdiag( y2 ))) + dMdm = sp.hstack((dMdm1, dMdm2, dMdm3, dMdm4, dMdm5, dMdm6)) + + return dMdm # ------------------------ Geometries ------------------------------ # diff --git a/SimPEG/Mesh/TensorMesh.py b/SimPEG/Mesh/TensorMesh.py index 5fd53344..e142d78d 100644 --- a/SimPEG/Mesh/TensorMesh.py +++ b/SimPEG/Mesh/TensorMesh.py @@ -549,25 +549,25 @@ class TensorMesh(BaseRectangularMesh, TensorView, DiffOperators, InnerProducts): return Utils.sdiag(Av.T * V * Utils.mkvc(materialProperty)) - def _fastFaceInnerProductDeriv(self, materialProperty=None): + def _fastFaceInnerProductDeriv(self, materialProperty=None, v=None): """ :param numpy.array materialProperty: 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) """ - return self._fastInnerProductDeriv('F', materialProperty=materialProperty) + return self._fastInnerProductDeriv('F', materialProperty=materialProperty, v=v) - def _fastEdgeInnerProductDeriv(self, materialProperty=None): + def _fastEdgeInnerProductDeriv(self, materialProperty=None, v=None): """ :param numpy.array materialProperty: 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', materialProperty=materialProperty) + return self._fastInnerProductDeriv('E', materialProperty=materialProperty, v=v) - def _fastInnerProductDeriv(self, AvType, materialProperty=None): + def _fastInnerProductDeriv(self, AvType, materialProperty=None, v=None): """ :param str AvType: 'E' or 'F' :param numpy.array materialProperty: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6)) @@ -576,11 +576,16 @@ class TensorMesh(BaseRectangularMesh, TensorView, DiffOperators, InnerProducts): """ if materialProperty is None or materialProperty.size == self.nC: Av = getattr(self, 'ave'+AvType+'2CC') - return self.dim * Av.T * Utils.sdiag(self.vol) + V = Utils.sdiag(self.vol) + if v is None: + return self.dim * Av.T * Utils.sdiag(self.vol) + return Utils.sdiag(v) * self.dim * Av.T * V if materialProperty.size == self.nC*self.dim: # anisotropic Av = getattr(self, 'ave'+AvType+'2CCV') V = sp.kron(sp.identity(self.dim), Utils.sdiag(self.vol)) - return Av.T * V + if v is None: + return Av.T * V + return Utils.sdiag(v) * Av.T * V if __name__ == '__main__': diff --git a/SimPEG/Tests/test_massMatricesDerivs.py b/SimPEG/Tests/test_massMatricesDerivs.py index 80f6bb83..8be5e4eb 100644 --- a/SimPEG/Tests/test_massMatricesDerivs.py +++ b/SimPEG/Tests/test_massMatricesDerivs.py @@ -5,45 +5,82 @@ from TestUtils import checkDerivative class TestInnerProductsDerivs(unittest.TestCase): - def setUp(self): - pass - def test_FaceIP_derivs_isotropic(self): - for d in range(3): - mesh = Mesh.TensorMesh([10,5,4][d:]) - v = np.random.rand(mesh.nF) - def fun(sig): - M = mesh.getFaceInnerProduct(sig) - Md = mesh.getFaceInnerProductDeriv(sig) - return M*v, Utils.sdiag(v)*Md - sig = np.random.rand(mesh.nC) - passed = checkDerivative(fun, sig, plotIt=False) - self.assertTrue(passed) + def doTestFace(self, h, rep, vec, fast): + mesh = Mesh.TensorMesh(h) + v = np.random.rand(mesh.nF) + 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(mesh.nC*rep) + return checkDerivative(fun, sig, num=5, plotIt=False) + + def doTestEdge(self, h, rep, vec, fast): + mesh = Mesh.TensorMesh(h) + v = np.random.rand(mesh.nE) + 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(mesh.nC*rep) + return checkDerivative(fun, sig, num=5, plotIt=False) + + def test_FaceIP_1D_isotropic(self): + self.assertTrue(self.doTestFace([10],1,True, False)) + def test_FaceIP_2D_isotropic(self): + self.assertTrue(self.doTestFace([10, 4],1,True, False)) + def test_FaceIP_3D_isotropic(self): + self.assertTrue(self.doTestFace([10, 4, 5],1,True, False)) + def test_FaceIP_2D_anisotropic(self): + self.assertTrue(self.doTestFace([10, 4],2,True, False)) + def test_FaceIP_3D_anisotropic(self): + self.assertTrue(self.doTestFace([10, 4, 5],3,True, False)) + def test_FaceIP_2D_tensor(self): + self.assertTrue(self.doTestFace([10, 4],3,True, False)) + def test_FaceIP_3D_tensor(self): + self.assertTrue(self.doTestFace([10, 4, 5],6,True, False)) + + def test_FaceIP_1D_isotropic_fast(self): + self.assertTrue(self.doTestFace([10],1, False, True)) + def test_FaceIP_2D_isotropic_fast(self): + self.assertTrue(self.doTestFace([10, 4],1, False, True)) + def test_FaceIP_3D_isotropic_fast(self): + self.assertTrue(self.doTestFace([10, 4, 5],1, False, True)) + def test_FaceIP_2D_anisotropic_fast(self): + self.assertTrue(self.doTestFace([10, 4],2, False, True)) + def test_FaceIP_3D_anisotropic_fast(self): + self.assertTrue(self.doTestFace([10, 4, 5],3, False, True)) - def test_EdgeIP_derivs_isotropic(self): - for h in [[10,5],[10,5,4]]: - mesh = Mesh.TensorMesh(h) - v = np.random.rand(mesh.nE) - def fun(sig): - M = mesh.getEdgeInnerProduct(sig) - Md = mesh.getEdgeInnerProductDeriv(sig) - return M*v, Utils.sdiag(v)*Md - sig = np.random.rand(mesh.nC) - passed = checkDerivative(fun, sig, plotIt=False) - self.assertTrue(passed) + def test_EdgeIP_2D_isotropic(self): + self.assertTrue(self.doTestEdge([10, 4],1,True, False)) + def test_EdgeIP_3D_isotropic(self): + self.assertTrue(self.doTestEdge([10, 4, 5],1,True, False)) + def test_EdgeIP_2D_anisotropic(self): + self.assertTrue(self.doTestEdge([10, 4],2,True, False)) + def test_EdgeIP_3D_anisotropic(self): + self.assertTrue(self.doTestEdge([10, 4, 5],3,True, False)) + def test_EdgeIP_2D_tensor(self): + self.assertTrue(self.doTestEdge([10, 4],3,True, False)) + def test_EdgeIP_3D_tensor(self): + self.assertTrue(self.doTestEdge([10, 4, 5],6,True, False)) + + def test_EdgeIP_2D_isotropic_fast(self): + self.assertTrue(self.doTestEdge([10, 4],1, False, True)) + def test_EdgeIP_3D_isotropic_fast(self): + self.assertTrue(self.doTestEdge([10, 4, 5],1, False, True)) + def test_EdgeIP_2D_anisotropic_fast(self): + self.assertTrue(self.doTestEdge([10, 4],2, False, True)) + def test_EdgeIP_3D_anisotropic_fast(self): + self.assertTrue(self.doTestEdge([10, 4, 5],3, False, True)) - def test_FaceIP_derivs_anisotropic(self): - for d in range(3): - mesh = Mesh.TensorMesh([10,5,4][d:]) - v = np.random.rand(mesh.nF) - def fun(sig): - M = mesh.getFaceInnerProduct(sig) - Md = mesh.getFaceInnerProductDeriv(sig) - return M*v, Utils.sdiag(v)*Md - sig = np.random.rand(mesh.nC*mesh.dim) - passed = checkDerivative(fun, sig, plotIt=False) - self.assertTrue(passed)