From 57f16fc408ba5829306839e7b5265581e1184861 Mon Sep 17 00:00:00 2001
From: rowanc1 <rowan@visiblegeology.com>
Date: Tue, 17 Jun 2014 14:01:57 -0600
Subject: [PATCH] updates (a bit slower implementation that is more general)

---
 simpegEM/TDEM/TDEM_b.py | 25 +++++++++++++++----------
 1 file changed, 15 insertions(+), 10 deletions(-)

diff --git a/simpegEM/TDEM/TDEM_b.py b/simpegEM/TDEM/TDEM_b.py
index ac1254f1..853746cb 100644
--- a/simpegEM/TDEM/TDEM_b.py
+++ b/simpegEM/TDEM/TDEM_b.py
@@ -107,13 +107,15 @@ class ProblemTDEM_b(BaseTDEMProblem):
         # fake initial 'e' fields
         p[:, 'e', 0] = 0.0
         tt = TensorType(self.mesh, self.curModel.transform)
-        c = self.mesh.getEdgeInnerProductDeriv(tt)()*(self.curModel.transformDeriv*vec)
+        dMdsig = self.mesh.getEdgeInnerProductDeriv(tt)
+        dsigdm_x_v = self.curModel.transformDeriv*vec
         for i in range(1,self.nT+1):
             # TODO: G[1] may be dependent on the model
             #       for a galvanic source (deriv of the dc problem)
             #
             # Do multiplication for all tx in self.survey.txList
-            p[:, 'e', i] = - sdiag(c) * u[:,'e',i] # i.e.: - sdiag(MsigDeriv * v) * e
+            for tx in self.survey.txList:
+                p[tx, 'e', i] = - dMdsig(u[tx,'e',i]) *  dsigdm_x_v
         return p
 
     def Gtvec(self, m, vec, u=None):
@@ -129,17 +131,20 @@ class ProblemTDEM_b(BaseTDEMProblem):
         if u is None:
             u = self.fields(m)
         self.curModel = m
+        tt = TensorType(self.mesh, self.curModel.transform)
+        dMdsig = self.mesh.getEdgeInnerProductDeriv(tt)
+        dsigdm = self.curModel.transformDeriv
 
-        nTx, nE = self.survey.nTx, self.mesh.nE
-        tmp = np.zeros(nE)
+        nTx = self.survey.nTx
+        VUs = None
         # Here we can do internal multiplications of Gt*v and then multiply by MsigDeriv.T in one go.
         for i in range(1,self.nT+1):
-            vu = vec[:,'e',i]*u[:,'e',i]
-            if nTx > 1:
-                vu = vu.sum(axis=1)
-            tmp += vu
-        tt = TensorType(self.mesh, self.curModel.transform)
-        p = -mkvc(self.curModel.transformDeriv.T*(self.mesh.getEdgeInnerProductDeriv(tt)().T*tmp))
+            vu = None
+            for tx in self.survey.txList:
+                vutx = dMdsig(u[tx,'e',i]).T * vec[tx,'e',i]
+                vu = vutx if vu is None else vu + vutx
+            VUs = vu if VUs is None else VUs + vu
+        p = -dsigdm.T*VUs
         return p
 
     def solveAh(self, m, p):