From 3eb8224e28609f0abc46d07a231e8895f17d2ffc Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Sat, 1 Mar 2014 06:54:30 -0800 Subject: [PATCH] Edge inner products are derivatives w/ testing. --- SimPEG/Mesh/InnerProducts.py | 30 +++++++++++-- SimPEG/Mesh/TensorMesh.py | 57 ++++++++++++++++++++++--- SimPEG/Tests/test_massMatricesDerivs.py | 15 ++++++- 3 files changed, 92 insertions(+), 10 deletions(-) diff --git a/SimPEG/Mesh/InnerProducts.py b/SimPEG/Mesh/InnerProducts.py index 034cdd01..266661a7 100644 --- a/SimPEG/Mesh/InnerProducts.py +++ b/SimPEG/Mesh/InnerProducts.py @@ -77,7 +77,6 @@ class InnerProducts(object): :rtype: scipy.csr_matrix :return: M, the inner product matrix (nF, nF) """ - fast = None if hasattr(self, '_fastFaceInnerProductDeriv'): @@ -97,8 +96,6 @@ class InnerProducts(object): P[0].T * sp.identity(n) * P[0] - - def getEdgeInnerProduct(self, materialProperty=None, returnP=False, invertProperty=False): """ :param numpy.array materialProperty: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6)) @@ -107,6 +104,14 @@ class InnerProducts(object): :rtype: scipy.csr_matrix :return: M, the inner product matrix (nE, nE) """ + fast = None + + if returnP is False and hasattr(self, '_fastEdgeInnerProduct'): + fast = self._fastEdgeInnerProduct(materialProperty=materialProperty, invertProperty=invertProperty) + + if fast is not None: + return fast + if invertProperty: materialProperty = invPropertyTensor(self, materialProperty) @@ -146,6 +151,25 @@ class InnerProducts(object): else: return A + + def getEdgeInnerProductDeriv(self, materialProperty=None, P=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) + """ + + fast = None + + if hasattr(self, '_fastEdgeInnerProductDeriv'): + fast = self._fastEdgeInnerProductDeriv(materialProperty=materialProperty) + + if fast is not None: + return fast + + raise NotImplementedError('Derivatives for the material property specified are not yet implemented.') + + # ------------------------ Geometries ------------------------------ # # diff --git a/SimPEG/Mesh/TensorMesh.py b/SimPEG/Mesh/TensorMesh.py index f55c6df6..5fd53344 100644 --- a/SimPEG/Mesh/TensorMesh.py +++ b/SimPEG/Mesh/TensorMesh.py @@ -503,35 +503,82 @@ class TensorMesh(BaseRectangularMesh, TensorView, DiffOperators, InnerProducts): :rtype: scipy.csr_matrix :return: M, the inner product matrix (nF, nF) """ + return self._fastInnerProduct('F', materialProperty=materialProperty, invertProperty=invertProperty) + + + def _fastEdgeInnerProduct(self, materialProperty=None, invertProperty=False): + """ + Fast version of getEdgeInnerProduct. + This does not handle the case of a full tensor materialProperty. + + :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 + :param bool invertProperty: inverts the material property + :rtype: scipy.csr_matrix + :return: M, the inner product matrix (nE, nE) + """ + return self._fastInnerProduct('E', materialProperty=materialProperty, invertProperty=invertProperty) + + + def _fastInnerProduct(self, AvType, materialProperty=None, invertProperty=False): + """ + Fast version of getFaceInnerProduct. + This does not handle the case of a full tensor materialProperty. + + :param numpy.array materialProperty: 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 invertProperty: inverts the material property + :rtype: scipy.csr_matrix + :return: M, the inner product matrix (nF, nF) + """ if materialProperty is None: materialProperty = np.ones(self.nC) if materialProperty.size == self.nC: if invertProperty: materialProperty = 1./materialProperty - Av = self.aveF2CC + Av = getattr(self, 'ave'+AvType+'2CC') V = Utils.sdiag(self.vol) return self.dim * Utils.sdiag(Av.T * V * materialProperty) if materialProperty.size == self.nC*self.dim: if invertProperty: materialProperty = 1./materialProperty - Av = self.aveF2CCV + 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)) - def _fastFaceInnerProductDeriv(self, materialProperty=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) + + + def _fastEdgeInnerProductDeriv(self, materialProperty=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) + + + def _fastInnerProductDeriv(self, AvType, materialProperty=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)) + :rtype: scipy.csr_matrix + :return: M, the inner product matrix (nF, nF) + """ if materialProperty is None or materialProperty.size == self.nC: - Av = self.aveF2CC + Av = getattr(self, 'ave'+AvType+'2CC') return self.dim * Av.T * Utils.sdiag(self.vol) if materialProperty.size == self.nC*self.dim: # anisotropic - Av = self.aveF2CCV + Av = getattr(self, 'ave'+AvType+'2CCV') V = sp.kron(sp.identity(self.dim), Utils.sdiag(self.vol)) return Av.T * V diff --git a/SimPEG/Tests/test_massMatricesDerivs.py b/SimPEG/Tests/test_massMatricesDerivs.py index f5484568..80f6bb83 100644 --- a/SimPEG/Tests/test_massMatricesDerivs.py +++ b/SimPEG/Tests/test_massMatricesDerivs.py @@ -11,7 +11,6 @@ class TestInnerProductsDerivs(unittest.TestCase): def test_FaceIP_derivs_isotropic(self): for d in range(3): mesh = Mesh.TensorMesh([10,5,4][d:]) - M,Ps = mesh.getFaceInnerProduct(returnP=True) v = np.random.rand(mesh.nF) def fun(sig): M = mesh.getFaceInnerProduct(sig) @@ -21,10 +20,22 @@ class TestInnerProductsDerivs(unittest.TestCase): passed = checkDerivative(fun, sig, plotIt=False) self.assertTrue(passed) + + 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_FaceIP_derivs_anisotropic(self): for d in range(3): mesh = Mesh.TensorMesh([10,5,4][d:]) - M,Ps = mesh.getFaceInnerProduct(returnP=True) v = np.random.rand(mesh.nF) def fun(sig): M = mesh.getFaceInnerProduct(sig)