From 880780fce0baa175c55fe18f41fa52208a86ba8a Mon Sep 17 00:00:00 2001
From: Lindsey Heagy <lindseyheagy@gmail.com>
Date: Tue, 4 Mar 2014 11:09:12 -0800
Subject: [PATCH] Jv implemented for single TX

---
 simpegEM/FDEM/DataFDEM.py                | 24 ++++++++++----
 simpegEM/FDEM/FDEM.py                    | 41 ++++++++++++++++--------
 simpegEM/Tests/test_forward_EMproblem.py |  6 ++--
 3 files changed, 48 insertions(+), 23 deletions(-)

diff --git a/simpegEM/FDEM/DataFDEM.py b/simpegEM/FDEM/DataFDEM.py
index bd9877b2..7343f819 100644
--- a/simpegEM/FDEM/DataFDEM.py
+++ b/simpegEM/FDEM/DataFDEM.py
@@ -1,4 +1,4 @@
-from SimPEG import Utils, np
+from SimPEG import Utils, np, sp
 from SimPEG.Data import BaseData
 from FieldsFDEM import FieldsFDEM
 
@@ -28,13 +28,23 @@ class DataFDEM(BaseData):
         BaseData.__init__(self, **kwargs)
         Utils.setKwargs(self, **kwargs)
 
-    def projectFields(self, u):
-        #TODO: this is hardcoded to 1Tx
-        return self.Qrx.dot(u.b[:,:,0].T).T
+    @property
+    def nRx(self):
+        return self.rxLoc.shape[0]
 
-    def projectFieldsAdjoint(self, d):
-        # TODO: fix this
-        pass
+    def projectFields(self, u):
+        P = sp.identity(self.prob.mesh.nE)
+        Pes = range(self.nFreq)
+        #TODO: this is hardcoded to 1Tx
+        for i, freqInd in enumerate(range(self.nFreq)):
+            e = u.get_e(freqInd)
+            Pes[i] = P*e
+        Pe = np.concatenate(Pes)
+        return Pe
+
+    def projectFieldsDeriv(self, u):
+        # TODO : more general
+        return sp.identity(self.prob.mesh.nE)
 
     ####################################################
     # Interpolation Matrices
diff --git a/simpegEM/FDEM/FDEM.py b/simpegEM/FDEM/FDEM.py
index 6580bac9..7436266d 100644
--- a/simpegEM/FDEM/FDEM.py
+++ b/simpegEM/FDEM/FDEM.py
@@ -1,5 +1,4 @@
-from SimPEG import Problem, Solver
-import numpy as np
+from SimPEG import Problem, Solver, np, sp
 from scipy.constants import mu_0
 from SimPEG.Utils import sdiag, mkvc
 from FieldsFDEM import FieldsFDEM
@@ -57,7 +56,7 @@ class ProblemFDEM_e(Problem.BaseProblem):
         # TODO: this will not work if tensor conductivity
         self._MeSigmaI = sdiag(1/self.MeSigma.diagonal())
         #TODO: assuming constant mu
-        self._MfMui = self.mesh.getFaceInnerProduct(1/mu_0)
+        self._MfMui = self.mesh.getFaceInnerProduct(1/mu_0*np.ones(self.mesh.nC))
 
     ####################################################
     # Internal Methods
@@ -65,7 +64,7 @@ class ProblemFDEM_e(Problem.BaseProblem):
 
     def getA(self, freqInd):
         """
-            :param int tInd: Time index
+            :param int fInd: Frequency index
             :rtype: scipy.sparse.csr_matrix
             :return: A
         """
@@ -99,9 +98,26 @@ class ProblemFDEM_e(Problem.BaseProblem):
 
 
     def Jvec(self, m, v, u=None):
+        # TODO: only 1 transmitter for now
+        # TODO: all P's the same
         if u is None:
-            u = self.fields(m)
-        raise NotImplementedError('Jvec todo!')
+           u = self.fields(m)
+
+        Jvs = range(self.data.nFreq)
+        P  = self.data.projectFieldsDeriv(u)
+
+        for i, freqInd in enumerate(range(self.data.nFreq)):
+            e = u.get_e(freqInd)
+            omega = self.data.omega[freqInd]
+            for txInd in self.data.nTx
+            b = 1j*omega*self.mesh.getEdgeInnerProductDeriv(m,v=e)*self.model.transformDeriv(m)*v
+            A = self.getA(freqInd)
+            Ab = Solver(A, options=self.solveOpts).solve(b)
+            Jvs[i] = -P*Ab
+
+        Jv = np.concatenate(Jvs)
+
+        return Jv
 
 
     def Jtvec(self, m, v, u=None):
@@ -110,8 +126,6 @@ class ProblemFDEM_e(Problem.BaseProblem):
         raise NotImplementedError('Jtvec todo!')
 
 
-
-
 if __name__ == '__main__':
     from SimPEG import *
     import simpegEM as EM
@@ -135,11 +149,11 @@ if __name__ == '__main__':
     model = Model.LogModel(mesh)
 
     opts = {'txLoc':0.,
-            'txType':'VMD_MVP',
-            'rxLoc': rxLoc,
-            'rxType':'bz',
-            'freq': np.logspace(0,3,4),
-            }
+        'txType':'VMD_MVP',
+        'rxLoc': rxLoc,
+        'rxType':'bz',
+        'freq': np.logspace(0,3,4),
+        }
     dat = EM.FDEM.DataFDEM(**opts)
 
     prb = EM.FDEM.ProblemFDEM_e(mesh, model)
@@ -160,3 +174,4 @@ if __name__ == '__main__':
 
 
 
+
diff --git a/simpegEM/Tests/test_forward_EMproblem.py b/simpegEM/Tests/test_forward_EMproblem.py
index e846e6c2..937f55da 100644
--- a/simpegEM/Tests/test_forward_EMproblem.py
+++ b/simpegEM/Tests/test_forward_EMproblem.py
@@ -17,8 +17,8 @@ class TDEM_bTests(unittest.TestCase):
         mesh = Mesh.Cyl1DMesh([hx,hy], -hy.sum()/2)
 
         active = mesh.vectorCCz<0.
-        model = Model.ActiveModel(mesh, active, -8, nC=mesh.nCz) 
-        model = Model.ComboModel(mesh, 
+        model = Model.ActiveModel(mesh, active, -8, nC=mesh.nCz)
+        model = Model.ComboModel(mesh,
                     [Model.LogModel, Model.Vertical1DModel, model])
 
         opts = {'txLoc':0.,
@@ -35,7 +35,7 @@ class TDEM_bTests(unittest.TestCase):
         self.sigma = np.ones(mesh.nCz)*1e-8
         self.sigma[mesh.vectorCCz<0] = 1e-1
         self.sigma = np.log(self.sigma[active])
-        
+
         self.prb.pair(self.dat)
 
     def test_analitic_b(self):