From bbe2a8563c45e1539cf0e279cffa83b2b202a235 Mon Sep 17 00:00:00 2001
From: rowanc1 <rowan@visiblegeology.com>
Date: Sat, 12 Jul 2014 11:43:45 -0500
Subject: [PATCH] change innerproducts to massMatrices

---
 simpegDC/BaseDC.py                       | 29 ++++++++++++++++++------
 simpegDC/Tests/test_forward_DCproblem.py |  9 ++++++++
 2 files changed, 31 insertions(+), 7 deletions(-)

diff --git a/simpegDC/BaseDC.py b/simpegDC/BaseDC.py
index d110f3b9..a4c41e60 100644
--- a/simpegDC/BaseDC.py
+++ b/simpegDC/BaseDC.py
@@ -95,7 +95,24 @@ class ProblemDC(Problem.BaseProblem):
         Utils.setKwargs(self, **kwargs)
 
 
-    deleteTheseOnModelUpdate = ['_A']
+    deleteTheseOnModelUpdate = ['_A', '_Msig', '_dMdsig']
+
+    @property
+    def Msig(self):
+        if getattr(self, '_Msig', None) is None:
+            sigma = self.curModel.transform
+            Av = self.mesh.aveF2CC
+            self._Msig = Utils.sdInv(Utils.sdiag(self.mesh.dim * Av.T * (1/sigma)))
+        return self._Msig
+
+    @property
+    def dMdsig(self):
+        if getattr(self, '_dMdsig', None) is None:
+            sigma = self.curModel.transform
+            Av = self.mesh.aveF2CC
+            dMdprop = self.mesh.dim * Utils.sdiag(self.Msig.diagonal()**2) * Av.T * Utils.sdiag(1./sigma**2)
+            self._dMdsig = lambda Gu: Utils.sdiag(Gu) * dMdprop
+        return self._dMdsig
 
     @property
     def A(self):
@@ -114,9 +131,7 @@ class ProblemDC(Problem.BaseProblem):
         if getattr(self, '_A', None) is None:
             D = self.mesh.faceDiv
             G = self.mesh.cellGrad
-            sigma = self.curModel.transform
-            Msig = self.mesh.getFaceInnerProduct(sigma, invProp=True, invMat=True)
-            self._A = D*Msig*G
+            self._A = D*self.Msig*G
             # Remove the null space from the matrix.
             self._A[-1,-1] /= self.mesh.vol[-1]
             self._A = self._A.tocsc()
@@ -168,7 +183,7 @@ class ProblemDC(Problem.BaseProblem):
         D = self.mesh.faceDiv
         G = self.mesh.cellGrad
         # Derivative of inner product, $\left(\mathbf{M}_{1/\sigma}^f\right)^{-1}$
-        dMdsig = self.mesh.getFaceInnerProductDeriv(sigma, invProp=True, invMat=True)
+        dMdsig = self.dMdsig
         # Derivative of model transform, $\deriv{\sigma}{\m}$
         dsigdm_x_v = self.curModel.transformDeriv * v
 
@@ -203,7 +218,7 @@ class ProblemDC(Problem.BaseProblem):
         D = self.mesh.faceDiv
         G = self.mesh.cellGrad
         A = self.A
-        Av_dm = self.mesh.getFaceInnerProductDeriv(sigma, invProp=True, invMat=True)
+        dMdsig = self.dMdsig
         mT_dm = self.mapping.deriv(m)
 
         dCdu = A.T
@@ -212,7 +227,7 @@ class ProblemDC(Problem.BaseProblem):
 
         Jtv = 0
         for i, ui in enumerate(u.T):  # loop over each column
-            Jtv += Av_dm( G * ui ).T * ( D.T * w[:,i] )
+            Jtv += dMdsig( G * ui ).T * ( D.T * w[:,i] )
 
         Jtv = - mT_dm.T * ( Jtv )
         return Jtv
diff --git a/simpegDC/Tests/test_forward_DCproblem.py b/simpegDC/Tests/test_forward_DCproblem.py
index cd4e2d9d..778b6a70 100644
--- a/simpegDC/Tests/test_forward_DCproblem.py
+++ b/simpegDC/Tests/test_forward_DCproblem.py
@@ -50,5 +50,14 @@ class DCProblemTests(unittest.TestCase):
         self.assertTrue(passed)
 
 
+    def test_massMatrices(self):
+        Gu = np.random.rand(self.mesh.nF)
+        def derChk(m):
+            self.p.curModel = m
+            return [self.p.Msig * Gu, self.p.dMdsig(Gu)]
+        passed = Tests.checkDerivative(derChk, self.m0, plotIt=False)
+        self.assertTrue(passed)
+
+
 if __name__ == '__main__':
     unittest.main()